*! Source of lmoremata.mlib


*! {smcl}
*! {marker mm_kern}{bf:mm_kern.mata}{asis}
*! version 1.0.3  18aug2020  Ben Jann
version 9.2

local kernels ///
 epanechnikov ///
 epan2 ///
 biweight ///
 triweight ///
 cosine ///
 gaussian ///
 parzen ///
 rectangle ///
 triangle

local klist tokens("`kernels'")

local kname `"(k!="" ? k : "epan2")"'

mata:

string scalar _mm_unabkern(|string scalar k)
{
	return(mm_strexpand(k, `klist', "epan2", 0, k+": invalid kernel"))
}

pointer scalar _mm_findkern(|string scalar k)
{
	pointer(real matrix function) scalar K
	K = findexternal("mm_kern_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kern_"+`kname'+"() not found")
	return(K)
}

pointer scalar _mm_findkint(|string scalar k)
{
	pointer(real matrix function) scalar K
	K = findexternal("mm_kint_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kint_"+`kname'+"() not found")
	return(K)
}

pointer scalar _mm_findkderiv(|string scalar k)
{
	pointer(real matrix function) scalar K
	K = findexternal("mm_kderiv_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kderiv_"+`kname'+"() not found")
	return(K)
}

pointer scalar _mm_findkdel0(|string scalar k)
{
	pointer(real scalar function) scalar K
	K = findexternal("mm_kdel0_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kdel0_"+`kname'+"() not found")
	return(K)
}

real matrix mm_kern(string scalar k, real matrix x)
{
	return((*_mm_findkern(_mm_unabkern(k)))(x))
}

real matrix mm_kint(string scalar k, real scalar r, |real matrix x)
{
	return(mm_callf(mm_callf_setup(_mm_findkint(_mm_unabkern(k)), args()-2, x), r))
}

real matrix mm_kderiv(string scalar k, real matrix x)
{
	return((*_mm_findkderiv(_mm_unabkern(k)))(x))
}

real matrix mm_kdel0(string scalar k)
{
	return((*_mm_findkdel0(_mm_unabkern(k)))())
}

// mm_kern_*: kernel function
// mm_kint_*(1,x): integral of K(z) from -infty to x
// mm_kint_*(2,x): integral of K(z)^2 from -infty to x
// mm_kint_*(2): kernel "roughness"
// mm_kint_*(3,x): integral of z*K(z) from -infty to x
// mm_kint_*(4,x): integral of z^2*K(z) from -infty to x
// mm_kint_*(4): kernel variance
// mm_kvar_*: variance
// mm_kderiv*: derivative of kernel
// mm_kdel0_*: canonical bandwidth

// (scaled) Epanechnikov kernel function
real matrix mm_kern_epanechnikov(real matrix x)
{
	return((.75:-.15*x:^2)/sqrt(5):*(abs(x):<sqrt(5)))
}
real matrix mm_kint_epanechnikov(real scalar r, |real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+(.75*x-.05*x:^3)/sqrt(5)):*(abs(x):<sqrt(5)) + (x:>=sqrt(5)))
	if (r==2 & args()<2) return(3/5^1.5)
	if (r==2)
	 return((.3/sqrt(5):+(.1125*x-.015*x:^3+.0009*x:^5)):*(abs(x):<sqrt(5)) + 3/5^1.5*(x:>=sqrt(5)))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.1875*sqrt(5):+.375/sqrt(5)*x:^2-.0375/sqrt(5)*x:^4):*(abs(x):<sqrt(5)))
	if (r==4 & args()<2) return(1)
	if (r==4)
	 return((.5:+.25/sqrt(5)*x:^3-.03/sqrt(5)*x:^5):*(abs(x):<sqrt(5)) + (x:>=sqrt(5)))
	else _error(3300)
}
//real scalar mm_kvar_epanechnikov() return(1)
real matrix mm_kderiv_epanechnikov(real matrix x)
{
	return(-0.3/sqrt(5)*x:*(abs(x):<sqrt(5)))
}
real scalar mm_kdel0_epanechnikov() return((3/5^1.5)^.2)


// Epanechnikov kernel function
real matrix mm_kern_epan2(real matrix x)
{
	return((.75:-.75*x:^2):*(abs(x):<1))
}
real matrix mm_kint_epan2(real scalar r, |real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+.75*x-.25*x:^3):*(abs(x):<1) + (x:>=1))
	if (r==2 & args()<2) return(.6)
	if (r==2)
	 return((.3:+.5625*x-.375*x:^3+.1125*x:^5):*(abs(x):<1) + .6*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.1875:+.375*x:^2-.1875*x:^4):*(abs(x):<1))
	if (r==4 & args()<2) return(.2) // 1/5
	if (r==4)
	 return((.1:+.25*x:^3-.15*x:^5):*(abs(x):<1) + .2*(x:>=1))
	else _error(3300)
}
//real scalar mm_kvar_epan2() return(.2)
real matrix mm_kderiv_epan2(real matrix x)
{
	return((-1.5*x):*(abs(x):<1))
}
real scalar mm_kdel0_epan2() return(15^.2)


// Biweight kernel function
real matrix mm_kern_biweight(real matrix x)
{
	return(.9375*(1:-x:^2):^2:*(abs(x):<1))
}
real matrix mm_kint_biweight(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+.9375*(x-2/3*x:^3+.2*x:^5)):*(abs(x):<1) + (x:>=1))
	if (r==2 & args()<2) return(5/7)
	if (r==2)
	 return((5/14:+225/256*x-75/64*x:^3+135/128*x:^5-225/448*x:^7+25/256*x:^9):*(abs(x):<1) + 5/7*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.15625:+.46875*x:^2-.46875*x:^4 +.15625*x:^6):*(abs(x):<1))
	if (r==4 & args()<2) return(1/7)
	if (r==4)
	 return((1/14:+.3125*x:^3-.375*x:^5+15/112*x:^7):*(abs(x):<1) + (x:>=1)/7)
	else _error(3300)
}
//real scalar mm_kvar_biweight() return(1/7)
real matrix mm_kderiv_biweight(real matrix x)
{
	return(3.75*(x:^3:-x):*(abs(x):<1))
}
real scalar mm_kdel0_biweight() return(35^.2)


// Triweight kernel function
real matrix mm_kern_triweight(real matrix x)
{
	return(1.09375*(1:-x:^2):^3:*(abs(x):<1))
}
real matrix mm_kint_triweight(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+1.09375*(x-x:^3+.6*x:^5-1/7*x:^7)):*(abs(x):<1) + (x:>=1))
	if (r==2 & args()<2) return(350/429)
	if (r==2)
	 return((175/429:+1225/1024*(x-2*x:^3+3*x:^5-20/7*x:^7+5/3*x:^9-6/11*x:^11+1/13*x:^13)):*(abs(x):<1) + 350/429*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-35/256:+35/64*x:^2-105/128*x:^4+35/64*x:^6-35/256*x:^8):*(abs(x):<1))
	if (r==4 & args()<2) return(1/9)
	if (r==4)
	 return((1/18:+35/96*x:^3-21/32*x:^5+15/32*x:^7-35/288*x:^9):*(abs(x):<1) + (x:>=1)/9)
	else _error(3300)
}
//real scalar mm_kvar_triweight() return(1/9)
real matrix mm_kderiv_triweight(real matrix x)
{
	return(-6.5625*x:*(1:-x:^2):^2:*(abs(x):<1))   // 210/32 = 6.5625
}
real scalar mm_kdel0_triweight() return((9450/143)^.2)


// Cosine trace
real matrix mm_kern_cosine(real matrix x)
{
	return((1:+cos(2*pi()*x)):*(abs(x):<.5))
}
real matrix mm_kint_cosine(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+x+sin(2*pi()*x)/(2*pi())):*(abs(x):<.5) + (x:>=.5))
	if (r==2 & args()<2) return(1.5)
	if (r==2)
	 return((.75:+1.5*x+sin(2*pi()*x)/pi()+cos(2*pi()*x):*sin(2*pi()*x)/(4*pi())):*(abs(x):<.5) + 1.5*(x:>=.5))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.125+.25/pi()^2:+cos(2*pi()*x)/(4*pi()^2)+x:*sin(2*pi()*x)/(2*pi())+.5*x:^2):*(abs(x):<.5))
	if (r==4 & args()<2) return(.5*(1/6-1/pi()^2))
	if (r==4)
	 return((1/24-.25/pi()^2:-sin(2*pi()*x)/(4*pi()^3)+x:*cos(2*pi()*x)/(2*pi()^2)+x:^2:*sin(2*pi()*x)/(2*pi())+x:^3/3):*(abs(x):<.5) + .5*(1/6-1/pi()^2)*(x:>=.5))
	else _error(3300)
}
//real scalar mm_kvar_cosine() return(.5*(1/6-1/pi()^2))
real matrix mm_kderiv_cosine(real matrix x)
{
	return(-2*pi()*sin(2*pi()*x):*(abs(x):<.5))
}
real scalar mm_kdel0_cosine() return((3/(.5*(1/6-1/pi()^2)^2))^.2)


// Gaussian kernel function
real matrix mm_kern_gaussian(real matrix x)
{
	return(normalden(x)) // alternatively: exp(-0.5*((x)^2))/sqrt(2*pi())
}
real matrix mm_kint_gaussian(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)            return(normal(x))
	if (r==2 & args()<2) return(.5/sqrt(pi())) // = normalden(0,sqrt(2))
	if (r==2)            return(.5/sqrt(pi()) * normal(sqrt(2)*x))
	if (r==3 & args()<2) return(0)
	if (r==3)            return(-1/sqrt(2*pi())*exp(-.5*x:^2))
	if (r==4 & args()<2) return(1)
	if (r==4)            return(-1/sqrt(2*pi())*x:*exp(-.5*x:^2) + normal(x))
	else _error(3300)
}
//real scalar mm_kvar_gaussian() return(1)
real matrix mm_kderiv_gaussian(real matrix x)
{
	return(-x:*normalden(x))
}
real scalar mm_kdel0_gaussian() return((1/(4*pi()))^.1)


// Parzen kernel function
real matrix mm_kern_parzen(real matrix x)
{
	return(((4/3:-8*x:^2+8*abs(x):^3):*(abs(x):<=.5)) +
	 ((8*(1:-abs(x)):^3/3):*(abs(x):>1/2:&abs(x):<=1)))
}
real matrix mm_kint_parzen(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((2/3:+8/3*(x+1.5*x:^2+x:^3+.25*x:^4)):*(x:>=-1:&x:<-.5) +
	  (.5:+4/3*x-8/3*x:^3-2*x:^4):*(x:>=-.5:&x:<0) +
	  (.5:+4/3*x-8/3*x:^3+2*x:^4):*(x:>=0:&x:<=.5) +
	  (1/3:+8/3*(x-1.5*x:^2+x:^3-.25*x:^4)):*(x:>.5:&x:<=1) + (x:>1))
	if (r==2 & args()<2) return(302/315)
	if (r==2)
	 return((64/63:+64/9*(x+3*x:^2+5*x:^3+5*x:^4+3*x:^5+x:^6+x:^7/7)):*(x:>=-1:&x:<-.5) +
	  (151/315:+16/9*(x-4*x:^3-3*x:^4+36/5*x:^5+12*x:^6+36/7*x:^7)):*(x:>=-.5:&x:<0) +
	  (151/315:+16/9*(x-4*x:^3+3*x:^4+36/5*x:^5-12*x:^6+36/7*x:^7)):*(x:>=0:&x:<=.5) +
	  (-2/35:+64/9*(x-3*x:^2+5*x:^3-5*x:^4+3*x:^5-x:^6+x:^7/7)):*(x:>.5:&x:<=1) + 302/315*(x:>1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-2/15:+4/3*x:^2+8/3*x:^3+2*x:^4+8/15*x:^5):*(x:>=-1:&x:<-.5) +
	  (-7/60:+2/3*x:^2-2*x:^4-1.6*x:^5):*(x:>=-.5:&x:<0) +
	  (-7/60:+2/3*x:^2-2*x:^4+1.6*x:^5):*(x:>=0:&x:<=.5) +
	  (-2/15:+4/3*x:^2-8/3*x:^3+2*x:^4-8/15*x:^5):*(x:>.5:&x:<=1))
	if (r==4 & args()<2) return(1/12)
	if (r==4)
	 return((2/45:+8/9*x:^3+2*x:^4+1.6*x:^5+4/9*x:^6):*(x:>=-1:&x:<-.5) +
	  (1/24:+4/9*x:^3-1.6*x:^5-4/3*x:^6):*(x:>=-.5:&x:<0) +
	  (1/24:+4/9*x:^3-1.6*x:^5+4/3*x:^6):*(x:>=0:&x:<=.5) +
	  (7/180:+8/9*x:^3-2*x:^4+1.6*x:^5-4/9*x:^6):*(x:>.5:&x:<=1) + (x:>1)/12)
	else _error(3300)
}
//real scalar mm_kvar_parzen() return(1/12)
real matrix mm_kderiv_parzen(real matrix x)
{
	return(((-16*x+24*sign(x):*abs(x):^2):*(abs(x):<=.5)) +
	 ((-8*sign(x):*(1:-abs(x)):^2):*(abs(x):>1/2:&abs(x):<=1)))
}
real scalar mm_kdel0_parzen() return((4832/35)^.2)


// Rectangle kernel function
real matrix mm_kern_rectangle(real matrix x)
{
	return(.5*(round(abs(x),1e-8):<1))
} // 1e-8 same as in kdensity.ado, version 2.3.6 26jun2000 (Stata 7)
real matrix mm_kint_rectangle(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+.5*x):*(round(abs(x),1e-8):<1) + (round(x,1e-8):>=1))
	if (r==2 & args()<2) return(.5)
	if (r==2)
	 return((.25:+.25*x):*(round(abs(x),1e-8):<1) + .5*(round(x,1e-8):>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.25:+.25*x:^2):*(round(abs(x),1e-8):<1))
	if (r==4 & args()<2) return(1/3)
	if (r==4)
	 return((1/6:+x:^3/6):*(round(abs(x),1e-8):<1) + (round(x,1e-8):>=1)/3)
	else _error(3300)
}
//real scalar mm_kvar_rectangle() return(1/3)
real matrix mm_kderiv_rectangle(real matrix x)
{
	return(J(rows(x), cols(x), 0)) // actually, the derivative is infty at -1 and -infty at 1
}
real scalar mm_kdel0_rectangle() return((9/2)^.2)


// Triangle kernel function
real matrix mm_kern_triangle(real matrix x)
{
	return((1:-abs(x)):*(abs(x):<1))
}
real matrix mm_kint_triangle(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+x+.5*x:^2):*(x:>-1:&x:<0) +
	  (.5:+x-.5*x:^2):*(x:>=0:&x:<1) + (x:>=1))
	if (r==2 & args()<2) return(2/3)
	if (r==2)
	 return((x+x:^2+(1:+x:^3)/3):*(x:>-1:&x:<0) +
	  (x-x:^2+(1:+x:^3)/3):*(x:>=0:&x:<1) + 2/3*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-1/6:+.5*x:^2+x:^3/3):*(x:>-1:&x:<0) +
	  (-1/6:+.5*x:^2-x:^3/3):*(x:>=0:&x:<1))
	if (r==4 & args()<2) return(1/6)
	if (r==4)
	 return((1/12:+x:^3/3+.25*x:^4):*(x:>-1:&x:<0) +
	  (1/12:+x:^3/3-.25*x:^4):*(x:>=0:&x:<1) + (x:>=1)/6)
	else _error(3300)
}
//real scalar mm_kvar_triangle() return(1/6)
real matrix mm_kderiv_triangle(real matrix x)
{
	return(-sign(x):*(abs(x):<1))
}
real scalar mm_kdel0_triangle() return(24^.2)

end


*! {smcl}
*! {marker mm_quantile}{bf:mm_quantile.mata}{asis}
*! version 2.0.3  01dec2020  Ben Jann
version 9.2
mata:

real matrix mm_quantile(real matrix X, | real colvector w,
    real matrix P, real scalar d, real scalar fw)
{
    real colvector o
    pointer (real matrix) scalar XX, ww
    
    if (args()<2) w = 1
    if (args()<3) P = (0, .25, .50, .75, 1)'
    if (args()<4) d = 2
    if (args()<5) fw = 0
    if (!anyof(0::9, d)) {
        display("{err}{it:def} must be an integer in [0,9]")
        _error(3300)
    }
    if (missing(X) | missing(w)) _error(3351)
    if (any(w:<0)) {
        display("{err}{it:w} must not be negative")
        _error(3300)
    }
    // drop zero frequency observations (i.e. observations for which w = 0) and
    // determine whether weights can be ignored
    XX = &X; ww = &w
    if (rows(w)!=1) {
        if (rows(w)!=rows(X)) _error(3200)
        if (anyof(w,0)) {
            o = select(1::rows(w), w)
            if (cols(o)==0) o = J(0,1,.) // select() may return 0x0
            XX = &(X[o,])
            ww = &(w[o])
        }
        // weights can be ignored if constant and d = 1 or 2
        if ((d==1 | d==2) & rows(*ww)) { // (*ww may be void)
            if (allof(*ww, (*ww)[1])) ww = &1
        }
        // weights can be ignored if w = 1 for all observations
        else if (allof(*ww,1)) ww = &1
    }
    else if (w==0) {
        XX = &(J(0,cols(X),.))
        ww = &(J(0,1,.))
    }
    // compute quantiles
    if (cols(X)==1 & cols(P)!=1 & rows(P)==1)
        return(_mm_quantile_sort(*XX, *ww, P', d, fw)')
    if (cols(X)!=1 & cols(P)!=1 & cols(X)!=cols(P)) _error(3200)
    return(_mm_quantile_sort(*XX, *ww, P, d, fw))
}

real matrix _mm_quantile_sort(real matrix X, real colvector w,
    real matrix P, real scalar d, real scalar fw)
{
    real scalar    i, c, c1, c2
    real colvector p, sX, sw, pP, sP
    real matrix    Q

    c1 = cols(X); c2 = cols(P)
    c = max((c1,c2))
    Q = J(rows(P), c, .)
    if (c1==c2) {
        if (w==1) {
            for (i=c; i; i--) {
                Q[,i] = _mm_quantile_d(sort(X[,i],1), editmissing(P[,i],1), d)
            }
            return(Q)
        }
        if (rows(w)==1) {
            for (i=c; i; i--) {
                pP = order(P[,i],1); sP = P[pP,i]
                Q[pP,i] = _mm_quantile_w(sort(X[,i],1), w, sP, d, fw)
            }
            return(Q)
        }
        for (i=c; i; i--) {
            p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]
            pP = order(P[,i],1); sP = P[pP,i]
            Q[pP,i] = _mm_quantile_w(sX, sw, sP, d, fw)
        }
        return(Q)
    }
    if (c1==1) {
        if (w==1) {
            sX = sort(X,1)
            for (i=c; i; i--) {
                Q[,i] = _mm_quantile_d(sX, editmissing(P[,i],1), d)
            }
            return(Q)
        }
        if (rows(w)==1) {
            sX = sort(X,1)
            for (i=c; i; i--) {
                pP = order(P[,i],1); sP = P[pP,i]
                Q[pP,i] = _mm_quantile_w(sX, w, sP, d, fw)
            }
            return(Q)
        }
        p = order((X,w),(1,2)); sX = X[p]; sw = w[p]
        for (i=c; i; i--) {
            pP = order(P[,i],1); sP = P[pP,i]
            Q[pP,i] = _mm_quantile_w(sX, sw, sP, d, fw)
        }
        return(Q)
    }
    if (c2==1) {
        if (w==1) {
            sP = editmissing(P,1)
            for (i=c; i; i--) {
                Q[,i] = _mm_quantile_d(sort(X[,i],1), sP, d)
            }
            return(Q)
        }
        if (rows(w)==1) {
            pP = order(P,1); sP = P[pP]
            for (i=c; i; i--) {
                Q[pP,i] = _mm_quantile_w(sort(X[,i],1), w, sP, d, fw)
            }
            return(Q)
        }
        pP = order(P,1); sP = P[pP]
        for (i=c; i; i--) {
            p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]
            Q[pP,i] = _mm_quantile_w(sX, sw, sP, d, fw)
        }
        return(Q)
    }
    _error(3200)
}

real matrix _mm_quantile(real colvector X, | real colvector w, 
    real matrix P, real scalar d, real scalar fw)
{   // X assumed sorted and non-missing
    // w assumed non-negative and non-missing
    real scalar    i
    real colvector o
    real matrix    Q
    pointer (real matrix) scalar XX, ww
    
    if (args()<2) w = 1
    if (args()<3) P = (0, .25, .50, .75, 1)'
    if (args()<4) d = 2
    if (args()<5) fw = 0
    // drop zero frequency observations (i.e. observations for which w = 0) and
    // determine whether weights can be ignored
    XX = &X; ww = &w
    if (rows(w)!=1) {
        if (rows(w)!=rows(X)) _error(3200)
        if (anyof(w,0)) {
            o = select(1::rows(w), w)
            if (cols(o)==0) o = J(0,1,.) // select() may return 0x0
            XX = &(X[o,])
            ww = &(w[o])
        }
        // weights can be ignored if constant and d < 3
        if (d<3 & rows(*ww)) { // (*ww may be void)
            if (allof(*ww, (*ww)[1])) ww = &1
        }
        // weights can be ignored if w = 1 for all observations
        else if (allof(*ww,1)) ww = &1
    }
    else if (w==0) {
        XX = &(J(0,1,.))
        ww = &(J(0,1,.))
    }
    // compute weighted quantiles: requires sorted p
    if (*ww!=1) {  // not rows(*ww)!=1 !
        if (rows(P)==1) {
            o = order(P',1)
            Q = o // just to dimension q
            Q[o] = _mm_quantile_w(*XX, *ww, P[o]', d, fw)
            return(Q')
        }
        Q = J(rows(P), cols(P), .)
        for (i=cols(P); i; i--) {
            o = order(P[,i],1)
            Q[o,i] = _mm_quantile_w(*XX, *ww, P[o,i], d, fw)
        }
        return(Q)
    }
    // compute unweighted quantiles
    if (rows(P)==1) return(_mm_quantile_d(*XX, editmissing(P',1), d)')
    Q = J(rows(P), cols(P), .)
    for (i=cols(P); i; i--) Q[,i] = _mm_quantile_d(*XX, editmissing(P[,i],1), d)
    return(Q)
}

real colvector _mm_quantile_d(real colvector X, real colvector p, real scalar d)
{   // X assumed sorted and non-missing
    // p assumed nonmissing
    real scalar     n, eps
    real colvector  pn, j, j1, h

    n = rows(X)
    if ((rows(p)*n)==0) return(J(rows(p), 1, .)) // no obs or rows(p)==0
    if (n==1)           return(J(rows(p), 1, X)) // only one obs
    if      (d==0) pn = p :* n
    else if (d==1) pn = p :* n
    else if (d==2) pn = p :* n
    else if (d==3) pn = p :* n :- .5
                                                // pn = a + p*(n + 1 - a - b)
    else if (d==4) pn = p :* n                  //      a = 0, b = 1
    else if (d==5) pn = p :* n :+ .5            //      a = b = 0
    else if (d==6) pn = p :* (n + 1)            //      a = b = 0 
    else if (d==7) pn = 1 :+ p :* (n - 1)       //      a = b = 1
    else if (d==8) pn = 1/3 :+ p :* (n + 1/3)   //      a = b = 1/3
    else if (d==9) pn = 3/8 :+ p :* (n + 1/4)   //      a = b = 3/8
    else {
        display("{err}{it:def} must be an integer in [0,9]")
        _error(3300)
    }
    if (d==0) return(X[mm_clip(floor(pn):+1, 1, n)])
    if (d==1) return(X[mm_clip(ceil(pn), 1, n)])
    if (d<=3) {
        j = floor(pn)
        if (d==2) h = ((pn:>j) :+ 1) / 2
        else      h = (pn:>j) :| mod(j,2)
    }
    else {
        eps = 4 * epsilon(1) // handle rounding error as in R's quantile()
        j = floor(pn :+ eps)
        h = pn - j
        h = h :* (abs(h):>=eps)
    }
    j1 = mm_clip(j:+1, 1, n)
    _mm_clip(j, 1, n)
    return((1:-h) :* X[j] :+ h :* X[j1])
}

real colvector _mm_quantile_w(real colvector X, real colvector w, real colvector p, 
    real scalar d, real scalar fw)
{   // X assumed sorted and non-missing
    // w assumed non-missing and *strictly* positive
    // p assumed sorted
    real scalar n
    real matrix W
    
    n = rows(X)
    if (n==0)       return(J(rows(p), 1, .))
    if (n==1)       return(J(rows(p), 1, X))
    if (rows(w)==0) return(J(rows(p), 1, .))
    // frequency weights or d<3 (for which type of weights is irrelevant)
    if (fw | d<3) {
        W = _mm_ecdf2(X, w, 0, 1) // W = (uniq X, runningsum(w))
        if (d==0) return(_mm_quantile_w_0(W[,1], W[,2], p))
        if (d==1) return(_mm_quantile_w_1(W[,1], W[,2], p))
        if (d==2) return(_mm_quantile_w_2(W[,1], W[,2], p))
        if (d==3) return(_mm_quantile_w_3(W[,1], W[,2], p))
        if (d==4) return(_mm_quantile_w_d(W[,1], W[,2], p,   0,   0))
        if (d==5) return(_mm_quantile_w_d(W[,1], W[,2], p,  .5,   0))
        if (d==6) return(_mm_quantile_w_d(W[,1], W[,2], p,   0,   1))
        if (d==7) return(_mm_quantile_w_d(W[,1], W[,2], p,   1,  -1))
        if (d==8) return(_mm_quantile_w_d(W[,1], W[,2], p, 1/3, 1/3))
        if (d==9) return(_mm_quantile_w_d(W[,1], W[,2], p, 3/8, 1/4))
        display("{err}{it:def} must be an integer in [0,9]")
        _error(3300)
    }
    // other cases
    if (rows(w)==1) W = (1::n) * w
    else {
        W = mm_colrunsum(w, 1, 1) // use quad precision
        if (W[n]>=.) W = J(n,1,.)
    }
    if (d==3) return(_mm_quantile_w_3b(X, W, p))
    else {
        if      (d==4) W =  W             /  W[n]
        else if (d==5) W = (W :- w/2)     /  W[n]
        else if (d==6) W =  W             / (W[n] + W[n]/n)
        else if (d==7) W = (W :- w)       / (W[n] - W[n]/n)
        else if (d==8) W = (W :- w/3)     / (W[n] + W[n]/(3*n))
        else if (d==9) W = (W :- (3/8)*w) / (W[n] + W[n]/(4*n))
        else {
            display("{err}{it:def} must be an integer in [0,9]")
            _error(3300)
        }
        return(mm_fastipolate(W, X, p, 1))
    }
}

real colvector _mm_quantile_w_0(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi
    real colvector P, q
    
    j = rows(W) 
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<=pi) break
        }
        q[i] = x[j]
    }
    return(q)
}

real colvector _mm_quantile_w_1(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi
    real colvector P, q
    
    j = rows(W) 
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        q[i] = x[j]
    }
    return(q)
}

real colvector _mm_quantile_w_2(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, r, pi
    real colvector P, q
    
    j = r = rows(W) 
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        if (W[j]==pi) {
            if (j==r) q[i] = x[j]
            else      q[i] = (x[j] + x[j+1])/2
        }
        else q[i] = x[j]
    }
    return(q)
}

real colvector _mm_quantile_w_3(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi, lo, up, d0, d1
    real colvector P, q
    
    j = rows(W)
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        // find j such that W[j-1] < pi <= W[j]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        // reached bottom
        if (j==1) { 
            q[|1\i|] = J(i,1, x[1])
            break
        }
        // if pi is larger than or equal to the next integer after W[j-1],
        // then use the upper x-value; decide between lower and upper x-value
        // only if pi is between W[j-1] and next integer
        lo = W[j-1]
        up = floor(lo) + 1
        if (pi>=up) q[i] = x[j]
        else {
            // set upper boundary to minimum between W[j] and next integer
            // after W[j-1]
            up = min( (W[j], up) )
            // obtain distances between pi and boundaries
            d0 = pi - lo
            d1 = up - pi
            // case 1: pi closer to upper boundary
            if (d1<d0) q[i] = x[j]
            // case 2: pi closer to lower boundary
            else if (d0<d1) q[i] = x[j-1]
            // case 3: equal distance
            else {
                // use x-value from lower boundary if lower boundary is integer 
                // and even; else use x-value from upper boundary; this is an
                // arbitrary decision
                if (!mod(lo,2)) q[i] = x[j-1]
                else            q[i] = x[j]
            }
        }
    }
    return(q)
}

real colvector _mm_quantile_w_3b(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi, d0, d1
    real colvector P, q
    
    j = rows(W)
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        // reached bottom
        if (j==1) { 
            q[|1\i|] = J(i,1, x[1])
            break
        }
        // obtain distances
        d0 = pi - W[j-1]
        d1 = W[j] - pi
        // case 1: pi closer to upper boundary
        if (d1<d0) q[i] = x[j]
        // case 2: pi closer to lower boundary
        else if (d0<d1) q[i] = x[j-1]
        // case 3: equal distance
        else {
            // use x-value from upper boundary if index of upper boundary is
            // even; else use x-value from lower boundary
            if (!mod(j,2)) q[i] = x[j]
            else           q[i] = x[j-1]
        }
    }
    return(q)
}

real colvector _mm_quantile_w_d(real colvector x, real colvector W, 
    real colvector p, real scalar a, real scalar b)
{
    real scalar    i, j, pi, lo, up, h, T
    real colvector q, pk
    
    j  = rows(W)
    T  = W[j] + b
    pk = (W :- a) / T
    i  = rows(p)
    q  = J(i,1,.)
    for (; i; i--) {
        pi = p[i]
        // find j such that pk[j-1] < pi <= pk[j]
        for (; j>1; j--) {
            if (pk[j-1]<pi) break
        }
        // reached bottom
        if (j==1) {
            q[|1\i|] = J(i,1, x[1])
            break
        }
        // if pi is larger than or equal to pk of next integer after W[j-1],
        // then use the upper x-value; interpolate only if pi is between 
        // pk[j-1] and pk of next integer
        up = (floor(W[j-1]) + 1 - a) / T
        if (pi>=up) q[i] = x[j]
        else {
            lo   = pk[j-1]
            up   = min( (pk[j], up) )
            h    = (pi - lo) / (up - lo)
            q[i] = (1-h) * x[j-1] + h * x[j]
        }
    }
    return(q)
}


end


*! {smcl}
*! {marker mm_median}{bf:mm_median.mata}{asis}
*! version 1.1.1  Ben Jann  03sep2020
version 9.2
mata:

real rowvector mm_median(real matrix X, | real colvector w)
{
    if (args()<2) w = 1
    return(mm_quantile(X, w, .5, 2))
}

real scalar _mm_median(real colvector X, | real colvector w)
{   // X assumed sorted and non-missing
    // w assumed non-missing and strictly positive
    if (args()<2) w = 1
    return(_mm_quantile(X, w, .5, 2))
}

end


*! {smcl}
*! {marker mm_iqrange}{bf:mm_iqrange.mata}{asis}
*! version 1.1.0  Ben Jann  13jul2020
version 9.2
mata:

real matrix mm_iqrange(real matrix X, | real colvector w, real scalar def, 
    real scalar fw)
{
    real matrix q

    if (args()<2) w = 1
    if (args()<3) def = 2
    if (args()<4) fw = 0
    q = mm_quantile(X, w, (.25 \ .75), def, fw)
    return(q[2,] - q[1,])
}

real scalar _mm_iqrange(real colvector X, | real colvector w, real scalar def, 
    real scalar fw)
{   // X assumed sorted and non-missing
    // w assumed non-missing and strictly positive
    real matrix q

    if (args()<2) w = 1
    if (args()<3) def = 2
    if (args()<4) fw = 0
    q = _mm_quantile(X, w, (.25 \ .75), def, fw)
    return(q[2] - q[1])
}

end


*! {smcl}
*! {marker mm_ecdf}{bf:mm_ecdf.mata}{asis}
*! version 1.1.0  10jul2020  Ben Jann
version 9.2
mata:

real matrix mm_ecdf(real matrix X, | real colvector w, real scalar mid, 
    real scalar nonorm, real scalar brk)
{
    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0
    if (args()<5) brk = 0
    return(mm_ranks(X, w, (brk!=0 ? 0 : 3), mid, (nonorm!=0 ? 0 : 1)))
}

real colvector _mm_ecdf(real colvector X, | real colvector w, real scalar mid, 
    real scalar nonorm, real scalar brk)
{
    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0
    if (args()<5) brk = 0
    return(_mm_ranks(X, w, (brk!=0 ? 0 : 3), mid, (nonorm!=0 ? 0 : 1)))
}

real matrix mm_ecdf2(real colvector X, | real colvector w, real scalar mid, 
    real scalar nonorm)
{
    real colvector p

    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0
    if (rows(X)==0) return(J(0,2,.))
    if (rows(w)==1) {
        p = order(X,1)
        return(_mm_ecdf2(X[p], w, mid, nonorm))
    }
    if (rows(w)!=rows(X)) _error(3200)
    p = order((X,w),(1,2))
    return(_mm_ecdf2(X[p], w[p], mid, nonorm))
}

real matrix _mm_ecdf2(real colvector X, | real colvector w, real scalar mid, 
    real scalar nonorm)
{
    real scalar    n
    real colvector p, cdf

    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0

    // compute running sum
    n = rows(X)
    if (rows(w)==1) cdf = (1::n) * w
    else {
        if (rows(w)!=n) _error(3200)
        // treat missing values as missing and use quad precision
        cdf = mm_colrunsum(w, 1, 1)
        if (cdf[n]>=.) cdf = J(n,1,.)
    }

    // remove ties (select last obs in each group)
    p = select(1::n, _mm_unique_tag(X, 1))
    cdf = cdf[p]
    
    // normalize
    if (!nonorm) cdf = cdf / cdf[rows(cdf)]
    
    // midrank adjustment
    if (mid) cdf = cdf :- mm_diff(0\cdf)/2
    
    // return result
    return(X[p], cdf)
}

end


*! {smcl}
*! {marker mm_relrank}{bf:mm_relrank.mata}{asis}
*! version 1.1.0  13jul2020  Ben Jann
version 9.2
mata:

real matrix mm_relrank(
    real matrix X,
    real colvector w,
    real matrix Y,
    | real scalar mid,
      real scalar nonorm, 
      real scalar brk,
      real colvector w2)
{
    if (args()<4) mid = 0
    if (args()<5) nonorm = 0
    if (args()<6) brk = 0
    if (args()<7) w2 = 1
    if (rows(w)!=1 & rows(w)!=rows(X)) _error(3200)
    if (rows(w2)!=1 & rows(w2)!=rows(Y)) _error(3200)
    if (cols(X)==1 & cols(Y)!=1 & rows(Y)==1 & rows(w2)==1)
        return(_mm_relrank_sort(X, w, Y', mid, nonorm, brk, w2)')
    return(_mm_relrank_sort(X, w, Y, mid, nonorm, brk, w2))
}

real matrix _mm_relrank_sort(
    real matrix X,
    real colvector w,
    real matrix Y,
    real scalar mid,
    real scalar nonorm, 
    real scalar brk,
    real colvector w2)
{
    real scalar    i, c, c1, c2
    real colvector p, sX, sw, p2, sY, sw2
    real matrix    R

    c1 = cols(X); c2 = cols(Y)
    c = max((c1,c2))
    R = J(rows(Y), c, .)
    if (c1==c2) {
        for (i=c; i; i--) {
            if (rows(w)==1) {; p = order(X[,i],1); sX = X[p,i]; sw = w; }
            else {; p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]; }
            if (rows(w2)==1) {; p2 = order(Y[,i],1); sY = Y[p2,i]; sw2 = w2; }
            else {; p2 = order((Y[,i],w2),(1,2)); sY = Y[p2,i]; sw2 = w2[p2]; }
            R[p2,i] = _mm_relrank(sX, sw, sY, mid, nonorm, brk, sw2)
        }
        return(R)
    }
    if (c1==1) {
        if (rows(w)==1) {; p = order(X,1); sX = X[p]; sw = w; }
        else {; p = order((X,w),(1,2)); sX = X[p]; sw = w[p]; }
        for (i=c; i; i--) {
            if (rows(w2)==1) {; p2 = order(Y[,i],1); sY = Y[p2,i]; sw2 = w2; }
            else {; p2 = order((Y[,i],w2),(1,2)); sY = Y[p2,i]; sw2 = w2[p2]; }
            R[p2,i] = _mm_relrank(sX, sw, sY, mid, nonorm, brk, sw2)
        }
        return(R)
    }
    if (c2==1) {
        if (rows(w2)==1) {; p2 = order(Y,1); sY = Y[p2]; sw2 = w2; }
        else {; p2 = order((Y,w2),(1,2)); sY = Y[p2]; sw2 = w2[p2]; }
        for (i=c; i; i--) {
            if (rows(w)==1) {; p = order(X[,i],1); sX = X[p,i]; sw = w; }
            else {; p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]; }
            R[p2,i] = _mm_relrank(sX, sw, sY, mid, nonorm, brk, sw2)
        }
        return(R)
    }
    _error(3200)
}

real colvector _mm_relrank(
    real colvector X,
    real colvector w,
    real colvector Y,
    | real scalar mid,
      real scalar nonorm, 
      real scalar brk,
      real colvector w2)
{
    real matrix cdf
    
    if (args()<4) mid = 0
    if (args()<5) nonorm = 0
    if (args()<6) brk = 0
    if (args()<7) w2 = 1
    if (rows(X)==0 | rows(Y)==0) return(J(rows(Y),1,.))
    
    // obtain CDF at unique values of X
    cdf = _mm_ecdf2(X, w, 0, nonorm)
    
    // compute relative ranks
    if (brk==0) {
        if (mid==0) return(_mm_relrank_1(Y, cdf[,1], cdf[,2]))
                    return(_mm_relrank_2(Y, cdf[,1], cdf[,2]))
    }
    if (rows(w2)==1) {
        if (mid==0) return(_mm_relrank_3(Y, cdf[,1], cdf[,2]))
                    return(_mm_relrank_4(Y, cdf[,1], cdf[,2]))
    }
    if (mid==0)     return(_mm_relrank_5(Y, w2, cdf[,1], cdf[,2]))
                    return(_mm_relrank_6(Y, w2, cdf[,1], cdf[,2]))
}

real colvector _mm_relrank_1(real colvector x, real colvector y, real colvector cdf)
{    // case 1: brk = 0, mid = 0
    real scalar    i, j, xi
    real colvector r
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) r[i] = cdf[j]
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_2(real colvector x, real colvector y, real colvector cdf)
{   // case 2: brk = 0, mid = 1
    real scalar    i, j, xi
    real colvector r, step
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            if (y[j]==xi) r[i] = cdf[j] - step[j]/2
            else          r[i] = cdf[j]
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_3(real colvector x, real colvector y, real colvector cdf)
{
    // case 3: brk = 1, mid = 0, no weights
    real scalar    i, j, k, xi
    real colvector r, step

    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            r[i] = cdf[j]
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) continue // no ties
                r[|k\i-1|] = cdf[j] :- step[j] :* (i:-(k::i-1)) / (i-k+1)
                i = k
            }
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_4(real colvector x, real colvector y, real colvector cdf)
{
    // case 4: brk = 1, mid = 1, no weights
    real scalar    i, j, k, xi
    real colvector r, step
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) {
                    r[i] = cdf[j] - step[j] * 0.5
                    continue
                }
                r[|k\i|] = cdf[j] :- step[j] :* ((i+.5):-(k::i)) / (i-k+1)
                i = k
            }
            else r[i] = cdf[j]
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_5(real colvector x, real colvector w, 
    real colvector y, real colvector cdf)
{
    // case 5: brk = 1, mid = 0, weighted
    real scalar    i, j, k, xi, W
    real colvector r, step, ww
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            r[i] = cdf[j]
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) continue // no ties
                ww = mm_colrunsum(w[|k\i|])
                W  = ww[rows(ww)]
                if (W==0) {
                    r[|k\i-1|] = cdf[j] :- step[j] :* (i:-(k::i-1)) / (i-k+1)
                }
                else {
                    ww = ww[|1\rows(ww)-1|]
                    r[|k\i-1|] = cdf[j] :- step[j] :* (W:-ww) / W
                }
                i = k
            }
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_6(real colvector x, real colvector w, 
    real colvector y, real colvector cdf)
{
    // case 6: brk = 1, mid = 1, weighted
    real scalar    i, j, k, xi, W
    real colvector r, step, ww
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) {
                    r[i] = cdf[j] - step[j] * 0.5
                    continue // no ties
                }
                ww = mm_colrunsum(w[|k\i|])
                W  = ww[rows(ww)]
                if (W==0) {
                    // if all observations have zero weight
                    r[|k\i|] = cdf[j] :- step[j] :* ((i+.5):-(k::i)) / (i-k+1)
                }
                else {
                    ww = ww - 0.5 * w[|k\i|]
                    r[|k\i|] = cdf[j] :- step[j] :* (W:-ww) / W
                }
                i = k
            }
            else r[i] = cdf[j]
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}

end


*! {smcl}
*! {marker mm_ranks}{bf:mm_ranks.mata}{asis}
*! version 1.1.0  19oct2020  Ben Jann
version 9.2
mata:

real matrix mm_ranks(real matrix X, | real colvector w,
 real scalar method, real scalar mid, real scalar norm)
{
    real scalar c, r, i
    real matrix result

    if (args()<2) w = 1
    if (args()<3) method = 0
    if (args()<4) mid = 0
    if (args()<5) norm = 0
    c = cols(X)
    r = rows(X)
    if (rows(w)!=1 & rows(w)!=r) _error(3200)
    if (c<1) return(J(r,0,.))
    if (c==1) return(_mm_ranks_sort(X, w, method, mid, norm))
    result = J(r, c, .)
    for (i=1; i<=c; i++) {
        result[,i] = _mm_ranks_sort(X[,i], w, method, mid, norm)
    }
    return(result)
}

real colvector _mm_ranks_sort(real colvector x, real colvector w,
 real scalar method, real scalar mid, real scalar norm)
{
    real scalar    I
    real colvector p, r

    I = rows(x)
    if (I==0) return(J(0,1,.))
    if      (method==4 & rows(w)!=1) p = order((x,w),(1,2))
    else if (method==0 | method==4)  p = order(x,1)
    else if (rows(w)!=1) p = order((x,w),(1,2)) // include w for stable results
    else p = order(x,1)
    r = __mm_ranks(x[p], (rows(w)!=1 ? w[p] : w), method, mid, norm)
    r[p] = r
    return(r)
}

real colvector _mm_ranks(real colvector x, | real colvector w,
 real scalar method, real scalar mid, real scalar norm) // sorted input assumed
{
    if (args()<2) w = 1
    if (args()<3) method = 0
    if (args()<4) mid = 0
    if (args()<5) norm = 0
    if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200)
    return(__mm_ranks(x, w, method, mid, norm))
}

real colvector __mm_ranks(real colvector x, real colvector w,
 real scalar method, real scalar mid, real scalar norm) // sorted input assumed
{
    real scalar    i, I, i0, i1
    real colvector ranks, R
    pointer scalar wR
    pragma unset   R

    // compute running sum
    I = rows(x)
    if (I==0) return(J(0,1,.))
    if (rows(w)!=1) {
        // treat missing values as missing and use quad precision
        ranks = mm_colrunsum(w, 1, 1)
        if (ranks[I]>=.) return(J(I,1,.))
    }
    else ranks = (1::I) * w
    
    // no special treatment of ties; use fast computations 
    if (method==0 | method==4) {
        if (mid & norm) return((ranks :- w/2) / ranks[I])
        if (norm)       return(ranks / ranks[I])
        if (mid)        return(ranks :- w/2)
        return(ranks)
    }
    
    // normalize
    if (norm) ranks = ranks / ranks[I]
    
    // handle ties
    if (method==3) {            // use lowest rank
        for(i=I-1; i>=1; i--) {
            if (x[i]==x[i+1]) ranks[i] = ranks[i+1]
        }
    }
    else if (method==1) {       // use highest rank
        for(i=2; i<=I; i++) {
            if (x[i]==x[i-1]) ranks[i] = ranks[i-1]
        }
    }
    else if (method==2) {       // use average rank
        i0 = i1 = 1
        if (rows(w)==1) wR = &1
        else wR = &R
        for(i=2; i<=I; i++) {
            if (x[i0]==x[i]) {
                i1++
                if (i<I) continue
            }
            if (i0<i1) {
                R = (i0 \ i1)
                ranks[|R|] = J(i1-i0+1,1,1) :* mean(ranks[|R|], w[|*wR|])
            }
            i0 = i1 = i
        }
    }
    else _error(3498,"ties = " + strofreal(method) + " not allowed")
    
    // apply midpoint adjustment
    if (mid) {
        i0 = i1 = 0
        for (i=1; i<=I; i++) {
            if (ranks[i]!=i0) {
                i1 = (i0 + ranks[i])/2
                i0 = ranks[i]
            }
            ranks[i] = i1
        }
    }
    
    // done
    return(ranks)
}

end


*! {smcl}
*! {marker mm_ddens}{bf:mm_ddens.mata}{asis}
*! version 1.0.1  12aug2020  Ben Jann
version 9.2
mata:

real matrix mm_ddens(real colvector X, | real colvector w, real vector minmax0,
    real scalar n0, real scalar h0, real scalar qui) // h0 will be replaced
{
    real scalar    n, N, h, lo, up
    real rowvector minmax
    real colvector AT, W, a
    
    // defaults
    if (args()<2) w  = 1
    if (n0>=.) n0 = 2^14
    if (h0<=0) _error(3300)
    if (args()<6) qui = 0
    
    // check input
    if (missing(X) | missing(w)) _error(3351)
    if (any(w:<0)) {
        display("{err}{it:w} must not be negative")
        _error(3300)
    }
    if (n0<2) _error(3300)
    
    // evaluation grid
    n = 2^ceil(ln(n0)/ln(2)) // round up to next power of 2
    if (length(minmax0)>=1) lo = minmax0[1]
    if (length(minmax0)>=2) up = minmax0[2]
    minmax = minmax(X)
    if (lo>=.) lo = minmax[1] - (minmax[2]-minmax[1])/10 // min - 10% of range
    if (up>=.) up = minmax[2] + (minmax[2]-minmax[1])/10 // min + 10% of range
    if (lo>minmax[1]) {
        display("{err}data smaller than lower bound not allowed")
        _error(3300)
    }
    if (up<minmax[2]) {
        display("{err}data larger than upper bound not allowed")
        _error(3300)
    }
    // bin data
    N  = mm_nobs(X, w)
    AT = rangen(lo, up, n)
    if (mm_issorted(X)) W = _mm_exactbin(X, w, rangen(lo, up, n+1)) / N
    else                W = mm_fastexactbin(X, w, rangen(lo, up, n+1)) / N
        // need to use exact binning because linear binning would introduce
        // some (non-vanishing) bias at the boundaries (doubling the first and
        // last grid count does not seem to help); a consequence of exact 
        // binning is that the density estimate will be slightly shifted/stretched
        // to the left; this error can be substantial if the grid size is small,
        // but it vanishes with increasing grid size
    // obtain discrete cosine transform of binned data
    a = Re( (1 \ 2 * exp(1i * (1::n-1) * pi() / (2*n)))
         :* fft(W[mm_seq(1,n-1,2)] \ W[mm_seq(n,2,2)]) )
    // compute bandwidth (unless provided by user)
    if (h0<.) h = (h0 / (up-lo))^2
    else {
        h = _mm_ddens_h(AT, W, N, (1::n-1):^2, (a[2::n]/2):^2, qui)
        h0 = sqrt(h) * (up-lo) // return optimal bandwidth
    }
    // smooth discrete cosine transform using bandwidth and back-transform
    a = a :* exp(-(0::n-1):^2 * pi()^2 * h/2)
    a = Re(invfft((n * exp(-1i * (0::n-1) * pi() / (2*n))) :* a)) / (up-lo)
    a[mm_seq(1, n, 2) \ mm_seq(2, n, 2)] = a[1::n/2 \ n::n/2+1] // reorder
    // return density estimate and grid
    return((a, AT))
}

real scalar _mm_ddens_h(real colvector AT, real colvector W, real scalar N, 
    real colvector I, real colvector a2, real scalar qui)
{
    real scalar s, hmin, h_os, ax, bx, rc, h
    
    // compute oversmoothed bandwidth
    s = _mm_iqrange(AT, W) / 1.349
    if (s<=0) s = sqrt(variance(AT, W*N))
    else      s = min((s, sqrt(variance(AT, W*N))))
    h_os = (s * (243/(35*N))^.2 * mm_kdel0_gaussian() / (AT[rows(AT)]-AT[1]))^2
    // minimum bandwidth given grid size
    hmin = (.5/(rows(AT)-1) * mm_kdel0_gaussian()/mm_kdel0_rectangle())^2
    // find optimal h
    bx = h_os
    ax = max((hmin, bx/2))
    rc = mm_root(h=., &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
    if (rc==2) { // move down
        bx = ax
        while (1) {
            if (bx<=hmin) break // cannot go below hmin
            ax = max((hmin, bx/2))
            rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
            if (rc==2) bx = ax // continue moving down
            else break // this also stops if there is a change in 
                      //  direction (rc==3) (in this case: h = current bx)
        }
    }
    else if (rc==3) { // move up
        ax = bx; bx = ax*1.5
        while (1) {
            rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
            if (rc==3) {; ax = bx; bx = ax*1.5; } // continue moving up
            else break // this also stops if there is a change in 
                       //  direction (rc==2) (in this case: h = current ax)
        }
    }
    if (h<=hmin) { // solution is smaller than hmin
        ax = h_os; bx = ax*1.5
        rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
        if (rc==2) h = . // moving up does not help
        else if (rc==3) { // continue moving up
            ax = bx; bx = ax*1.5
            while (1) {
                rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
                if (rc==3) {; ax = bx; bx = ax*1.5; } // continue moving up
                else break // this also stops if there is a change in 
                           //  direction (rc==2) (in this case: h = current ax)
            }
        }
    }
    // return result
    if (h>=.) {
        if (qui==0) display("{txt}(bandwidth estimation failed; using oversmoothed rule)")
        return(h_os)
    }
    return(h)
}

real scalar _mm_ddens_fixpnt(real scalar h, real scalar N, real colvector I, 
    real colvector a2)
{
    real scalar    l, s, K0, c
    real colvector f, t
    
    l = 7
    f = 2 * pi()^(2*l) * sum(I:^l :* a2 :* exp(-I * pi()^2 * h))
    for (s=l-1; s>=2; s--) {
        K0 = mm_prod(mm_seq(1, 2*s-1, 2)) / sqrt(2*pi())
        c  = (1 + (1/2)^(s + 1/2)) / 3
        t  = (2 * c * K0/N/f):^(2/(3 + 2*s))
        f  = 2 * pi()^(2*s) * sum(I:^s :* a2 :* exp(-I * pi()^2 * t))
    }
    return((2 * N * sqrt(pi()) * f)^(-2/5) - h)
}

end



*! {smcl}
*! {marker mm_freq}{bf:mm_freq.mata}{asis}
*! version 1.0.4, Ben Jann, 09oct2020
version 9.1
mata:

real colvector mm_freq(transmorphic matrix x,
 | real colvector w, transmorphic matrix levels)
{
	real colvector p

	if (args()<2) w = 1
	if (args()<3) levels = .
	if (cols(x)==0) return(_mm_freq(x, w, levels))
	if (rows(w)==1) return(_mm_freq(sort(x,1..cols(x)), w, levels))
	p = order(x,1..cols(x))
	return(_mm_freq(x[p,], w[p,], levels))
}

real colvector _mm_freq(transmorphic matrix x,
 | real colvector w, transmorphic matrix levels)
{
	real scalar    i, j, l
	real colvector result

	if (args()<2) w = 1
	if (args()<3) levels = .
	if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200)
	if (levels==.) levels = _mm_uniqrows(x)
	if (rows(x)==0) return(J(0,1, .))
	l = rows(levels)
	result = J(l,1,0)
	j = 1
	for (i=1; i<=rows(x); i++) {
		for (;j<=l;j++) {
			if (x[i,]<=levels[j,]) break
		}
		if (j>l) break
		if (x[i,]!=levels[j,]) continue
		result[j] = result[j] + (rows(w)!=1 ? w[i] : w)
	}
	return(result)
}

real colvector mm_freq2(transmorphic matrix x,
 | real colvector w)
{
    real colvector p

    if (args()<2) w = 1
    if (cols(x)==0) return(_mm_freq2(x, w))
    p = order(x,1..cols(x))
    if (rows(w)==1) return(_mm_freq2(x[p,],w)[invorder(p)])
    return(_mm_freq2(x[p,],w[p,])[invorder(p)])
}

real colvector _mm_freq2(transmorphic matrix x,
 | real colvector w)
{
    real scalar    i, j
    real colvector result

    if (args()<2) w = 1
    if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200)
    if (rows(x)==0) return(J(0, 1, .))
    result = J(rows(x),1,.)
    j = 1
    for (i=2; i<=rows(x); i++) {
        if (x[i,]!=x[i-1,]) {
            result[|j \ i-1|] = J(i-j, 1, (rows(w)==1 ?
                (i-j)*w : sum(w[|j \ i-1|])))
            j = i
        }
    }
    result[|j \ i-1|] = J(i-j, 1, (rows(w)==1 ?
        (i-j)*w : sum(w[|j \ i-1|])))
    return(result)
}

end


*! {smcl}
*! {marker mm_histogram}{bf:mm_histogram.mata}{asis}
*! version 1.0.3, Ben Jann, 22jun2006
version 9.0
mata:

real matrix mm_histogram(
 real colvector x,     // data
 | real colvector w,  // weights
   real colvector g,  // grid (interval borders) (default: 10 intervals)
   real scalar dir)   // 0 right inclusive (default), else left inclusive
{
	real matrix H
	real scalar i

	if (args()<2) w = 1
	if (args()<3) g = rangen(min(x), max(x), 11)
	if (args()<4) dir = 0

	H = J(rows(g)-1, 3, 0)
	for (i=1; i<=rows(H); i++) {
		H[i,2] = g[i+1]-g[i]     // width of bins
		H[i,1] = g[i] + H[i,2]/2 // bin midpoints
	}
	H[.,3] = mm_exactbin(x, w, g, dir) :/ (mm_nobs(x,w) * H[.,2])
	return(H)
}
end


*! {smcl}
*! {marker mm_collapse}{bf:mm_collapse.mata}{asis}
*! version 1.0.1  18feb2009  Ben Jann
version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real matrix mm_collapse(
    real matrix    X,              // data
    real colvector w,              // weights (or 1)
    real colvector ID,             // subgroup IDs
    |  pointer(function) scalar f, // function; default: mean()
      `opts')                      // arguments to pass through to f()
{
    real colvector p
    pointer scalar Xs, IDs, ws

    if ((rows(X)!=rows(ID)) | (rows(w)!=1 & rows(X)!=rows(w))) _error(3200)
    p     = order(ID, 1)
    if (isfleeting(X))      Xs  = &(X = X[p,])
    else                    Xs  = &(X[p,])
    if (rows(w)==1)         ws  = &w
    else if (isfleeting(w)) ws  = &(w = w[p])
    else                    ws  = &(w[p])
    if (isfleeting(ID))     IDs = &(ID = ID[p])
    else                    IDs = &(ID[p])

    if      (args()==3)  return(_mm_collapse(*Xs, *ws, *IDs))
    else if (args()==4)  return(_mm_collapse(*Xs, *ws, *IDs, f))
    else if (args()==5)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1))
    else if (args()==6)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2))
    else if (args()==7)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3))
    else if (args()==8)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4))
    else if (args()==9)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5))
    else if (args()==10) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6))
    else if (args()==11) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7))
    else if (args()==12) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7, o8))
    else if (args()==13) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7, o8, o9))
    else if (args()==14) return(_mm_collapse(*Xs, *ws, *IDs, f, `opts'))
}


real matrix _mm_collapse(
    real matrix    X,               // data
    real colvector w,               // weights (or 1)
    real colvector ID,              // subgroup IDs
    |  pointer(function) scalar f0, // function; default: mean()
      `opts')                       // arguments to pass through to f()
{
    real scalar                 i, j, c, a, b, ww
    real matrix                 info, res, key
    pointer(function) scalar    f
    transmorphic                setup

    if (rows(X)!=rows(ID))     _error(3200)
    ww      = rows(w)!=1
    if (ww & rows(w)!=rows(X)) _error(3200)
    if (rows(X)<1)             return(J(0, cols(X)+1, missingof(X)))

    f       = args()<4 ? &mean() : f0
    setup = mm_callf_setup(f, args()-4, `opts')

    info = _mm_panels(ID)
    key  = J(rows(info), 1, missingof(ID))
    res  = J(rows(info), c=cols(X), missingof(X))
    b = 0
    for (i=1; i<=rows(info); i++) {
        a = b + 1
        b = b + info[i]
        key[i] = ID[a]
        for (j=1; j<=c; j++) {
            res[i, j] =
                mm_callf(setup, X[|a,j \ b,j|], ww ? w[|a \ b|] : w)
        }
    }
    return(key, res)
}

end


*! {smcl}
*! {marker mm_gini}{bf:mm_gini.mata}{asis}
*! version 1.0.3  26mar2008  Ben Jann
version 9.2
mata:

real rowvector mm_gini(real matrix X, | real colvector w)
{
    real rowvector result
    real scalar c, i

    if (args()==1) w = 1
    c = cols(X)
    if (c<1) return(J(1,0,.))
    if (c==1) return(_mm_gini(X, w))
    result = J(1, c, .)
    for (i=1; i<=c; i++) {
        result[1, i] = _mm_gini(X[,i], w)
    }
    return(result)
}

real scalar _mm_gini(real colvector x, real colvector w)
{
    real matrix mv

    if (rows(x)<1) return(.)
    mv = mm_meanvariance0((x, mm_ranks(x, w, 3, 1, 1)), w)
    return(mv[3,1] * 2 / mv[1,1])
}

end


*! {smcl}
*! {marker mm_nobs}{bf:mm_nobs.mata}{asis}
*! version 1.0.0, Ben Jann, 07jun2006
version 9.2
mata:

real scalar mm_nobs(transmorphic matrix x, real colvector w)
{
	return(rows(w)==1 ? rows(x)*w :
	 (rows(x)==rows(w) ? colsum(w) : _error(3200)))
}
end


*! {smcl}
*! {marker mm_greedy}{bf:mm_greedy.mata}{asis}
*! version 1.0.1  17may2019  Ben Jann
version 9.2
local Int   real scalar
local IntC  real colvector
local IntM  real matrix
local RS    real scalar
local RC    real colvector
local RM    real matrix
local T     transmorphic
local TM    transmorphic matrix
local PidC  pointer(`Int') colvector
local PidM  pointer(`Int') matrix
local dist  pointer(function) scalar
mata:

`RM' mm_greedy2(`TM' T, `TM' C, `Int' n, `RS' calip, `dist' f, | `T' fopts)
{
    return(mm_greedy_pairs(mm_greedy(T, C, n, calip, f, fopts)))
}

`RM' mm_greedy_pairs(`IntM' P)
{
    `Int'  i, j, k, n
    `IntC' N
    `RM'   E
    
    N = rownonmissing(P)
    k = sum(N)
    E = J(k, 3, .)
    for (i=rows(P); i; i--) {
        n = N[i]
        for (j=n; j; j--) {
            E[k,] = (i, P[i,j], 1/n)
            k--
        }
    }
    return(E)
}

`IntM' mm_greedy(`TM' T, `TM' C, `Int' n, `RS' calip, `dist' f, | `T' fopts)
{
    `RC'   d
    `PidM' ij
    pragma unset d
    
    // check input
    if (cols(T)!=cols(C)) _error(3200, "T and C must have same number of columns")
    if (cols(T)==0) _error(3200, "T and C must have at least one column")
    // make ID lookup table and compute distances
    if (calip>=.) ij = _mm_greedy_dist(d, T, C, f, fopts)
    else          ij = _mm_greedy_dist_calip(d, T, C, f, fopts, calip)
    // match and return IDs of matched controls
    if (n<=1 | n>=.) return(_mm_greedy_match_1(ij, d, rows(T)))
    return(_mm_greedy_match_n(ij, d, rows(T), n))
}

`PidM' _mm_greedy_dist(`RC' d, `TM' T, `TM' C, `dist' f, `T' fopts)
{
    `Int'  i, a, b, nT, nC, N
    `PidC' Ti, Ci
    `PidM' ij
    
    nT = rows(T); nC = rows(C)   // number of treated and controls
    N = nT * nC                  // number of combinations
    Ti = _mm_greedy_pid(nT)      // make ID pointers for treated
    Ci = _mm_greedy_pid(nC)      // make ID pointers for controls
    d = J(N, 1, .)               // prepare vector of distances
    ij = J(N, 2, NULL)           // prepare ID lookup table
    if (N==0) return(ij)         // T or C is void
    b = 0
    for (i=nT; i; i--) {         // loop over treated => compare to all controls
        a = b + 1                // start index of current batch
        b = b + nC               // end index of current batch
        d[|a\b|] = (*f)(T[i,], C, fopts)   // add distances to distance vector
        ij[|a,1\b,.|] = J(nC,1,Ti[i]), Ci  // add ID pointers to lookup table
    }
    return(ij)
}

`PidM' _mm_greedy_dist_calip(`RC' d, `TM' T, `TM' C, `dist' f, `T' fopts,
    `RS' calip)
{
    `Int'  i, a, b, nT, nC, N, n
    `PidC' Ti, Ci
    `PidM' ij
    `RC'   D
    `IntC' P, p
    
    nT = rows(T); nC = rows(C)    // number of treated and controls
    N = nT * nC                   // number of combinations
    Ti = _mm_greedy_pid(nT)       // make ID pointers for treated
    Ci = _mm_greedy_pid(nC)       // make ID pointers for controls
    d = J(N, 1, .)                // prepare vector of distances
    ij = J(N, 2, NULL)            // prepare ID lookup table
    if (N==0) return(ij)          // T or C is void
    b = 0; P = (1::nC)
    for (i=nT; i; i--) {          // loop over treated => compare to all controls
        D = (*f)(T[i,], C, fopts) // compute distances
        p = select(P, D:<=calip)  // index of distances within caliper
        n = length(p)             // number of distances within caliper
        if (n==0) continue        // no valid distances -> skip case
        a = b + 1                 // start index of current batch
        b = b + n                 // end index of current batch
        if (n==nC) {              // all distances are valid
            d[|a\b|] = D          // add distances to distance vector
            ij[|a,1\b,.|] = J(nC,1,Ti[i]), Ci // add ID pointers to lookup table
            continue
        }
        d[|a\b|] = D[p]          // add distances to distance vector
        ij[|a,1\b,.|] = J(n,1,Ti[i]), Ci[p] // add ID pointers to lookup table
    }
    if (b==N) return(ij)         // all distances valid; all rows filled
    if (b==0) {                  // no valid distances
        d = J(0, 1, .)
        return(J(0, 2, NULL))
    }
    d = d[|1\b|]                 // remove unused rows from d
    return(ij[|1,1\b,.|])        // remove unused rows from ij
}

`PidC' _mm_greedy_pid(`Int' n)
{
    `Int'  i
    `PidC' P
    
    P = J(n, 1, NULL)                 // prepare container for ID pointers
    for (i=n; i; i--) P[i] = &(i*1)   // generate ID pointers
    return(P)
}

`IntC' _mm_greedy_match_1(`PidM' ij, `RC' d, `Int' nT)
{
    `Int'  i, j, ii, jj, N, L
    `IntC' p, M
    
    M = J(nT, 1, .)                        // prepare container for matched IDs
    L = nT                                 // maximum number of matches
    N = rows(d)                            // number of distances
    p = order(d, 1)                        // obtain sort order of distances
    for (i=1; i<=N; i++) {
        j = p[i]                           // get next position
        if ((ii = *(ij[j,1]))==.) continue // treatment case already matched
        if ((jj = *(ij[j,2]))==.) continue // control case already used
        M[ii] = jj                         // copy matched control ID
        *(ij[j,1]) = .                     // delete treatment ID
        *(ij[j,2]) = .                     // delete control ID
        if (!(--L)) break                  // maximum reached; no treatment IDs left
    }
    return(M)
}

`IntM' _mm_greedy_match_n(`PidM' ij, `RC' d, `Int' nT, `Int' n)
{
    `Int'  i, j, ii, jj, N, L, nc
    `IntC' p, NC
    `IntM' M
    
    NC = J(nT, 1, 0)                       // counter for number of matches
    M = J(nT, n, .)                        // container for matched IDs
    L = nT * n                             // maximum number of matches
    N = rows(d)                            // number of distances
    p = order(d, 1)                        // obtain sort order of distances
    for (i=1; i<=N; i++) {
        j = p[i]                           // get next position
        if ((ii = *(ij[j,1]))==.) continue // treatment case already matched
        if ((jj = *(ij[j,2]))==.) continue // control case already used
        nc = NC[ii]; nc++                  // update number of matches
        M[ii, nc] = jj                     // copy matched control ID
        if (nc==n) *(ij[j,1]) = .          // delete treatment ID
        *(ij[j,2]) = .                     // delete control ID
        NC[ii] = nc                        // store number of matches
        if (!(--L)) break                  // maximum reached; no treatment IDs left
    }
    return(M)
}

end




*! {smcl}
*! {marker mm_colvar}{bf:mm_colvar.mata}{asis}
*! version 1.0.1, Ben Jann, 22jun2006
version 9.0
mata:

real rowvector mm_colvar(real matrix X, |real colvector w)
{
	real rowvector result
	real scalar j, k

	if (args()==1) w = 1

	result = J(1, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[j] = variance(X[,j], w)
	}
	return(result)
}

real matrix mm_meancolvar(real matrix X, |real colvector w)
{
	real matrix result
	real scalar j, k

	if (args()==1) w = 1

	result = J(2, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[.,j] = meanvariance(X[,j], w)
	}
	return(result)
}

end


*! {smcl}
*! {marker mm_variance0}{bf:mm_variance0.mata}{asis}
*! version 1.0.1, Ben Jann, 22jun2006
version 9.0
mata:

real matrix mm_variance0(real matrix X, |real colvector w)
{
	real rowvector	CP
	real rowvector	means
	real scalar	n

	if (args()==1) w = 1

	CP = quadcross(w,0, X,1)
	n  = cols(CP)
	means = CP[|1\n-1|] :/ CP[n]
	if (missing(CP) | CP[n]==0) return(J(cols(X), cols(X), .))
	return(crossdev(X,0,means, w, X,0,means) :/ CP[n])
}

real matrix mm_meanvariance0(real matrix X, |real colvector w)
{
	real rowvector	CP
	real rowvector	means
	real scalar	n

	if (args()==1) w = 1

	CP = quadcross(w,0, X,1)
	n  = cols(CP)
	means = CP[|1\n-1|] :/ CP[n]
	if (missing(means)) return(means \ J(cols(X),cols(X),.))
	return(means \ crossdev(X,0,means, w, X,0,means) :/ CP[n])
}

end


*! {smcl}
*! {marker mm_mse}{bf:mm_mse.mata}{asis}
*! version 1.0.0, Ben Jann, 22jun2006
version 9.0
mata:

real matrix mm_mse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector CP

	CP = cross(w,0, X,1)
	if (missing(CP) | CP[cols(CP)]==0) return(J(cols(X), cols(X), .))
	return(mm_sse(X, w, mu) :/ CP[cols(CP)])
}

real matrix mm_sse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector CP

	CP = cross(w,0, X,1)
	if (missing(CP) | CP[cols(CP)]==0) return(J(cols(X), cols(X), .))
	return(crossdev(X,0,mu, w, X,0,mu))
}

real rowvector mm_colmse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector result
	real scalar j, k

	result = J(1, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[j] = mm_mse(X[,j], w, mu[j])
	}
	return(result)
}

real rowvector mm_colsse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector result
	real scalar j, k

	result = J(1, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[j] = mm_sse(X[,j], w, mu[j])
	}
	return(result)
}

end


*! {smcl}
*! {marker mm_sample}{bf:mm_sample.mata}{asis}
*! version 1.0.2, Ben Jann, 23oct2020
version 9.1
mata:

real colvector mm_sample(
 real colvector n,
 real matrix strata,
 | real colvector cluster,
   real colvector w,
   real scalar wor,
   real scalar count,
   real scalar fast,
   real scalar alt,
   real scalar nowarn)
{
	real scalar i, b, e, n0, n1, N, cb, ce, ups
	real colvector s, si, ci, wi, nn, R

	if (args()<3) cluster=.
	if (args()<4) w=1
	if (args()<5) wor=0
	if (args()<6) count=0
	if (args()<7) fast=0
	if (args()<8) alt=0
	if (args()<9) nowarn=0
	if (cols(strata)<1) _error(3200)
	ups = (rows(w)!=1)
	if (fast==0) {
		if (rows(strata)>1 | (cluster==. & ups==0)) {
			if (missing(strata)) _error(3351, "'strata' has missing values")
		}
	}
	if (rows(n)<1) _error(3200)

//check w
	if (ups) {
		if (cluster!=.) {
			if (rows(w)!=rows(cluster))
				_error(3200, "rows('w') unequal number of clusters")
		}
		if (fast==0) {
			if (missing(w)) _error(3351, "'w' has missing values")
			if (colsum(w:<0)) _error(3498, "'w' has negative values")
			if (colsum(w)<=0) _error(3498, "sum('w') is zero")
		}
	}

//simple sample, unstratified
	if (rows(strata)==1) {
		N = strata[1,1]
		if (cluster==.) {
			if (N>=.) N = rows(w)
			else if (ups) {
				if (N!=rows(w))
					_error(3200, "rows('w') unequal population size")
			}
			return(_mm_sample(n, (ups ? w : N), ups, wor, count, alt, nowarn))
		}

//cluster sample, unstratified
		if (rows(cluster)<1) _error(3200)
		if (N>=.) N = colsum(cluster)
		else if (fast==0) {
			if (N!=colsum(cluster))
			  _error(3200, "sum of cluster sizes unequal population size")
		}
		s = _mm_sample(n, (ups ? w : rows(cluster)), ups, wor, count, alt, nowarn)
		return(_mm_expandclusters(s, cluster, count, N))
	}

//stratified sample
	if (rows(strata)<1) _error(3200)
	if (cluster!=. & cols(strata)<2) _error(3200)
//-sample sizes
	if (rows(n)==1) {
		if (n>=.) {
			if (cluster==.) nn = strata[.,1]
			else            nn = strata[.,2]
		}
		else nn = J(rows(strata),1,n)
	}
	else {
		if (rows(n)!=rows(strata)) _error(3200)
		nn = n
		for (i=1;i<=rows(nn);i++) {
			if (nn[i]>=.) {
				if (cluster==.) nn[i] = strata[i,1]
				else            nn[i] = strata[i,2]
			}
		}
	}
//-prepare sample vector
	if (count) s = J(colsum(strata[.,1]),1,.)
	else s = J(colsum(nn),1,.)
//-draw samples within strata
	if (count==0) n0 = 1
	e = 0
	ce = 0
	for (i=1;i<=rows(strata);i++) {
		b = e + 1
		e = e + strata[i,1]
//--simple sample
		if (cluster==.) {
			si = _mm_sample(nn[i], (ups ? w[|b \ e|] : strata[i,1]), ups, wor,
				count, alt, nowarn)
		}
//--cluster sample
		else {
			cb = ce + 1
			ce = ce + strata[i,2]
			ci = cluster[|cb \ ce|]
			wi = (ups ? w[|cb \ ce|] : rows(ci))
			if (fast==0) {
				if (strata[i,1]!=colsum(ci))
				  _error(3200, "sum of cluster sizes unequal size of stratum")
			}
			si = _mm_sample(nn[i], wi, ups, wor, count, alt, nowarn)
			if (count) si = _mm_expandclusters(si, ci, 1, strata[i,1])
		}
//--add subsample to sample vector
		if (count) R = (b \ e)
		else {
			n1 = n0 + rows(si) - 1
			R = (n0 \ n1)
			if (cluster==.) si = si :+ (b-1)
			else si = si :+ (cb-1)
			n0 = n1 + 1
		}
		if (R[1]<=R[2]) s[|R|] = si
	}
	if (count==0&cluster!=.) s = _mm_expandclusters(s, cluster)
	return(s)
}

real colvector _mm_sample(real scalar n, real colvector NorW, 
    real scalar ups, real scalar wor, real scalar count, real scalar alt, 
    real scalar nowarn)
{
    if (ups) {
        if (wor) return(mm_upswor(n, NorW, count, nowarn))
        return(mm_upswr(n, NorW, count))
    }
    if (wor) return(mm_srswor(n, NorW, count, alt))
    return(mm_srswr(n, NorW, count))
} 

real colvector _mm_expandclusters(real colvector s0,
 real colvector cluster, | real scalar count, real scalar N)
{
	real scalar i, e, b, n0, n1
	real colvector s, eclust

	if (args()<3) count = 0
	if (count) {
		s = J(N,1,.)
		e = 0
		for (i=1;i<=rows(cluster);i++) {
			b = e + 1
			e = e + cluster[i]
			s[|b \ e|] = J(e-b+1, 1, s0[i])
		}
		return(s)
	}
	s = J(colsum(cluster[s0,]),1,.)
	eclust = mm_colrunsum(cluster)
	n0 = 1
	for (i=1;i<=rows(s0);i++) {
		e = eclust[s0[i]]
		b = e - cluster[s0[i]] + 1
		n1 = n0 + e-b
		s[|n0 \ n1|] = (b::e)
		n0 = n1+1
	}
	return(s)
}

end


*! {smcl}
*! {marker mm_srswr}{bf:mm_srswr.mata}{asis}
*! version 1.0.0, Ben Jann, 27mar2006
version 9.1
mata:

// simple random sample, with replacement
real colvector mm_srswr(real scalar n, real scalar N,
 | real scalar count)
{
	real colvector res, u
	real scalar i, nn

// check args
	if (args()<3) count = 0
	if (N>=.) _error(3351)
	if (N<1) _error(3300)
	if (n>=.) nn = N
	else nn = n

// no sampling
	if (N==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}

// draw sample
	u = ceil(uniform(nn,1)*N)
	if (count==0) return(u)
	res = J(N,1,0)
	for (i=1;i<=nn;i++) {
		res[u[i]] = res[u[i]] + 1
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_srswor}{bf:mm_srswor.mata}{asis}
*! version 1.0.1, Ben Jann, 23oct2020
version 9.1
mata:

// simple random sample, without replacement
real colvector mm_srswor(real scalar n, real scalar N,
 | real scalar count, real scalar alt)
{
	real colvector res, u
	real scalar i, nn

// check args
	if (args()<3) count = 0
	if (args()<4) alt = 0
	if (N>=.) _error(3351)
	if (N<1) _error(3300)
	if (n<0) _error(3300)
	if (n>=.) nn = N
	else nn = n
	if (N<nn) _error(3300, "n may not be larger than N")

// no sampling
	if (N==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}
	if (nn==0) {
		if (count) return(J(N,1,0))
		return(J(0,1,.))
	}

// draw sample
	if (alt) {
		if (nn^2 < N/2) u = _mm_srswor_a(n, N)
		else            u = _mm_srswor_b(n, N)
	}
	else u = mm_unorder2(N)[|1 \ nn|] //=> stable results even for very large N
	if (count==0) return(u)
	res = J(N,1,0)
	for (i=1;i<=nn;i++) {
		res[u[i]] = 1
	}
	return(res)
}

real colvector _mm_srswor_a(real scalar n, real scalar N)
{   // naive algorithm; fast when n is small compared to N
    real scalar    i, u
    real colvector H
    
    H = J(n, 1, .)
    H[1] = ceil(uniform(1,1)*N)
    for (i=2; i<=n; i++) {
        while (anyof(H[|1 \ i|], u = ceil(uniform(1,1)*N))) continue
        H[i] = u
    }
    return(H)
}

real colvector _mm_srswor_b(real scalar n, real scalar N)
{   // Fisher–Yates shuffle; fast when n is large compared to N
    real scalar    i, j, k
    real colvector I, H
    
    k = N
    I = 1::k
    H = J(n, 1, .)
    for (i=1; i<=n; i++) {
        j = ceil(uniform(1,1)*k)
        H[i] = I[j]
        I[j] = I[k]
        k--
    }
    return(H)
}


end


*! {smcl}
*! {marker mm_upswr}{bf:mm_upswr.mata}{asis}
*! version 1.0.0, Ben Jann, 01apr2006
version 9.1
mata:

// unequal probability sampling, with replacement
real colvector mm_upswr(real scalar n, real colvector w,
 | real scalar count)
{
	real colvector ub, res, u
	real scalar i, j, nn

// check args
	if (args()<3) count = 0
	if (rows(w)<1) _error(3498, "no cases")
	if (n>=.) nn = rows(w)
	else nn = n

// no sampling
	if (rows(w)==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}

// draw sample
	ub = mm_colrunsum(w)
	ub = ub/ub[rows(ub)]
	u = sort(uniform(nn,1),1)
	j=1
	if (count) {
		res = J(rows(w),1,0)
		for (i=1;i<=nn;i++) {
			while (u[i]>ub[j]) j++
			res[j] = res[j]+1
		}
		return(res)
	}
	res = J(nn,1,0)
	for (i=1;i<=nn;i++) {
		while (u[i]>ub[j]) j++
		res[i] = j
	}
	return(mm_jumble2(res)) //=> stable results even for very large n
}

end


*! {smcl}
*! {marker mm_upswor}{bf:mm_upswor.mata}{asis}
*! version 1.0.1, Ben Jann, 12jun2007
version 9.1
mata:

// unequal probability sampling, without replacement
real colvector mm_upswor(real scalar n, real colvector w,
 | real scalar count, real scalar nowarn)
{
	real colvector p, ub, res
	real scalar i, j, nn, u, z

// check args
	if (args()<4) nowarn = 0
	if (args()<3) count = 0
	if (rows(w)<1) _error(3498, "no cases")
	if (n<0) _error(3300)
	if (n>=.) nn = rows(w)
	else nn = n
	if (nowarn==0) {
		if (rows(w)<nn) _error(3300, "n may not be larger than rows(w)")
	}

// no sampling
	if (rows(w)==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}

// draw sample
	p = mm_unorder2(rows(w)) //=> stable results even for very large N
	ub = w[p] / colsum(w) * nn
	if (nowarn==0) {
		z = colsum(ub:>1)
		if (z==1) _error(3300, "1 case has w_i*n/sum(w)>1")
		if (z>0) _error(3300, strofreal(z) + " cases have w_i*n/sum(w)>1")
	}
	ub = mm_colrunsum(ub)
	u = uniform(1,1)
	j=1
	if (count) {
		res = J(rows(w),1,0)
		for (i=1;i<=nn;i++) {
			while (u>ub[j]) j++
			res[j] = res[j]+1
			u = u + 1
		}
		return(res[invorder(p)])
	}
	res = J(nn,1,0)
	for (i=1;i<=nn;i++) {
		while (u>ub[j]) j++
		res[i] = p[j]
		u = u + 1
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_jk}{bf:mm_jk.mata}{asis}
*! version 1.0.0, Ben Jann, 04jul2006
version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

struct mm_jkstats {
	real colvector reps, ndrop, fpc
	real rowvector stat
	real matrix    rstat
}

// note: reps is total n. of replications
//       ndrop is n. of unsuccessful replications
//       rows(rstat) = reps-ndrop is n. of successful replications

struct mm_jkstats scalar mm_jk(
 pointer(real rowvector function) scalar f,
 real matrix X,
 | real colvector w,
   real scalar nodots,       // 0
   real colvector strata,    // .
   real colvector cluster,   // .
   real colvector subpop,    // .
   real colvector fpc,       // .
   real matrix stat,         // .
   `opts')                   // opts to pass to f
{
	real scalar               i, ii, i0, i1, h, h0, h1, j
	real colvector            C, wjk, wjki
	real matrix               S
	transmorphic              fsetup
	struct mm_jkstats scalar  jk
	pointer scalar            ws

	if (args()<3) w = 1
	if (args()<4) nodots = 0
	if (args()<7) subpop = .
	if (args()<8) fpc = .
	if (args()<9) stat = .

	if (subpop!=. & rows(subpop)>0) ws = &(w :* (subpop:!=0))
	else ws = &w
	if (rows(*ws)==1) ws = &J(rows(X), 1, *ws)

	fsetup = mm_callf_setup(f, args()-9, `opts')

	if (stat==.) jk.stat = mm_callf(fsetup, X, *ws)[1,]
	else         jk.stat = stat[1,]

	mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.)
	if (anyof(S, 1)) _error(460, "singleton cluster detected")

	if (nodots==0) {
		printf("{txt}\nJackknife replications ({res}%g{txt})\n",
		 colsum(S[,2]))
		display("{txt}{hline 4}{c +}{hline 3} 1 " +
		 "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " +
		 "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ")
	}
	jk.rstat = J(colsum(S[,2]), cols(jk.stat), .)
	jk.ndrop = J(rows(S), 1, 0)
	jk.fpc   = J(rows(S), 1, 0)
	h1 = i1 = 0
	i = ii = 1
	for (h=1; h<=rows(S); h++) { // for each stratum
		h0 = h1 + 1
		h1 = h1 + S[h,1]
		if (fpc!=. & rows(fpc)>0) jk.fpc[h] = fpc[h1]
		wjk = *ws
		wjk[| h0 \ h1 |] = wjk[| h0 \ h1 |] * (S[h,2] / (S[h,2] - 1))
		for (j=1; j<=S[h,2]; j++) { // for each PSU
			i0 = i1 + 1
			i1 = i1 + (C!=. ? C[i] : 1)
			wjki = wjk
			wjki[| i0 \ i1 |] = J(i1-i0+1,1,0)
			jk.rstat[ii,] = mm_callf(fsetup, X, wjki)[1,]
			if (missing(jk.rstat[ii,])) {
				jk.ndrop[h] = jk.ndrop[h] + 1
				if (nodots==0) printf("{err}x{txt}")
			}
			else {
				ii++
				if (nodots==0) printf(".")
			}
			if (nodots==0) {
				if (mod(i,50)==0) printf(" %5.0f\n",i)
				displayflush()
			}
			i++
		}
	}
	if (nodots==0 & mod(i-1,50)) display("")
	if (ii<i) {
		if (ii>1) jk.rstat = jk.rstat[|1,1 \ ii-1,.|]
		else jk.rstat = J(0, cols(jk.rstat), .)
	}
	jk.reps = S[,2]
	return(jk)
}

real matrix mm_jk_report(
 struct mm_jkstats scalar jk,    // jackknife estimates
 | string vector what,           // b,pseudo,mean,bias,v,se,ci
   real scalar level,            //
   real scalar mse,              // !=0 => use mse formula
   real vector fpc0)             // sampling fraction (for each stratum)
{
	real scalar    theta, pseudo, mean, bias, se, v, ci, l,
	               i, alpha, tval
	real rowvector jkse, jkv, jkmean
	real colvector fpc, nh, mh
	real matrix    res, jkpseudo

	if (args()<2) what  = "se"
	if (args()<3) level = st_numscalar("c(level)")
	if (args()<4) mse   = 0

	theta = pseudo = mean = bias = v = se = ci = l = 0
	for (i=1; i<=length(what); i++) {
		if      (what[i]=="theta" | what[i]=="b") theta = (l = l + 1)
		else if (what[i]=="pseudo") pseudo = (l = l + rows(jk.rstat))
		else if (what[i]=="mean"  ) mean  = (l = l + 1)
		else if (what[i]=="bias"  ) bias  = (l = l + 1)
		else if (what[i]=="v"     ) v     = (l = l + cols(jk.rstat))
		else if (what[i]=="se"    ) se    = (l = l + 1)
		else if (what[i]=="ci"    ) ci    = (l = l + 2)
		else _error(3498, "'" + what[i] + "' not allowed")
	}

	res = J(l, cols(jk.rstat), .)
	nh = jk.reps - jk.ndrop
	fpc = fpc0
	if (cols(fpc)!=1) fpc = fpc'
	if (rows(fpc)==0 | fpc==.) fpc = jk.fpc
	if (rows(fpc)==0)          fpc = J(rows(jk.reps),1,0)
// observed parameter vector
	if (theta) res[theta,] = jk.stat
// pseudo values
	if (pseudo | (mse==0 & (mean | bias))) {
		jkpseudo = jk.rstat +
		 mm_expand(nh, nh, 1, 1):*(jk.stat :- jk.rstat)
		if (pseudo)
		 res[|pseudo-rows(jk.rstat)+1,1 \ pseudo,.|] = jkpseudo
	}
// jackknife mean
	if (mean | bias) {
		if (mse) jkmean = mean(jk.rstat, 1)
		else     jkmean = mean(jkpseudo, 1)
		if (mean) res[mean,] = jkmean
	}
// bias
	if (bias) res[bias,] = jkmean - jk.stat
// variance
	if (v) {
		mh = mm_expand((1:-fpc) :* (nh :- 1) :/ nh, nh, 1, 1)
		if (mse) jkv = mm_sse(jk.rstat, mh, jk.stat)
		else     jkv = mm_sse(jk.rstat, mh, mean(jk.rstat,1))
		res[|v-cols(res)+1,1 \ v,.|] = jkv
	}
// standard errors
	if (se | ci) {
		if (v) jkse = sqrt(diagonal(jkv))'
		else {
			mh = mm_expand((1:-fpc) :* (nh :- 1) :/ nh, nh, 1, 1)
			if (mse) jkse = sqrt(mm_colsse(jk.rstat, mh, jk.stat))
			else     jkse = sqrt(mm_colsse(jk.rstat, mh, mean(jk.rstat,1)))
		}
		if (se) res[se,] = jkse
	}
// confidence interval
	if (ci) {
		alpha = (100-level)/200
		tval  = invttail(rows(jk.rstat)-rows(jk.reps), alpha)
		res[ci-1,] = jk.stat - tval*jkse
		res[ci,]   = jk.stat + tval*jkse
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_bs}{bf:mm_bs.mata}{asis}
*! version 1.0.0, Ben Jann, 10jul2006
version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

struct mm_bsstats {
	real scalar    reps, ndrop
	real rowvector stat, se
	real matrix    rstat, rse
}

// note: reps is total n. of replications
//       ndrop is n. of unsuccessful replications
//       rows(rstat) = reps-ndrop is n. of successful replications

struct mm_bsstats scalar mm_bs(
 pointer(real rowvector function) scalar f,
 real matrix X,
 | real colvector w,
   real scalar reps,        // 50
   real scalar d,           // 0
   real scalar nodots,      // 0
   real colvector strata,   // .
   real colvector cluster,  // .
   real matrix stat,        // .
   `opts')                  // opts to pass to f
{
	real scalar               i, ii, hasse
	real colvector            C, N, s
	real matrix               S, res
	transmorphic              setup
	struct mm_bsstats scalar  bs

	if (args()<3) w = 1
	if (args()<4) reps = 50
	if (args()<5) nodots = 0
	if (args()<6) d = 0
	if (args()<9) stat = .

	setup = mm_callf_setup(f, args()-9, `opts')

	if (stat==.) {
		res = mm_callf(setup, X, w)
		bs.stat = res[1,]
		if (hasse=(rows(res)>1)) bs.se = res[2,]
	}
	else {
		bs.stat = stat[1,]
		if (hasse=(rows(stat)>1)) bs.se = stat[2,]
	}

	mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.)
	if (anyof(S, 1)) _error(460, "singleton cluster detected")
	if (d==0|d>=.) N = .
	else N = S[,2] :- d

	if (nodots==0) {
		printf("{txt}\nBootstrap replications ({res}%g{txt})\n", reps)
		display("{txt}{hline 4}{c +}{hline 3} 1 " +
		 "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " +
		 "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ")
	}
	bs.rstat = J(reps, cols(bs.stat), .)
	if (hasse) bs.rse = J(reps, cols(bs.se), .)
	for (i=ii=1; i<=reps; i++) {
		s  = mm_sample(N, S, C)
		res = mm_callf(setup, X[s,], rows(w)==1 ? w : w[s,])
		if (missing(res)) {
			if (nodots==0) printf("{err}x{txt}")
		}
		else {
			bs.rstat[ii,] = res[1,]
			if (hasse) bs.rse[ii,] = res[2,]
			ii++
			if (nodots==0) printf(".")
		}
		if (nodots==0) {
			if (mod(i,50)==0) printf(" %5.0f\n",i)
			displayflush()
		}
	}
	if (nodots==0 & mod(i-1,50)) display("")
	if (ii<i) {
		if (ii>1) {
			bs.rstat = bs.rstat[|1,1 \ ii-1,.|]
			if (hasse) bs.rse = bs.rse[|1,1 \ ii-1,.|]
		}
		else {
			bs.rstat = J(0, cols(bs.rstat), .)
			if (hasse) bs.rse = J(0, cols(bs.rse), .)
		}
	}
	bs.reps = reps
	bs.ndrop = reps - rows(bs.rstat)
	return(bs)
}

struct mm_bsstats scalar mm_bs2(
 pointer(real rowvector function) scalar f,
 real matrix X,
 | real colvector w,
   real scalar reps,        // 50
   real scalar d,           // 0
   real scalar nodots,      // 0
   real colvector strata,   // .
   real colvector cluster,  // .
   real matrix stat,        // .
   `opts')                  // opts to pass to f
{
	real scalar               i, ii, hasse
	real colvector            C, N, ws
	real matrix               S, res
	transmorphic              setup
	struct mm_bsstats scalar  bs

	if (args()<3) w = 1
	if (args()<4) reps = 50
	if (args()<5) nodots = 0
	if (args()<6) d = 0
	if (args()<9) stat = .

	setup = mm_callf_setup(f, args()-9, `opts')

	if (stat==.) {
		res = mm_callf(setup, X, w)
		bs.stat = res[1,]
		if (hasse=(rows(res)>1)) bs.se = res[2,]
	}
	else {
		bs.stat = stat[1,]
		if (hasse=(rows(stat)>1)) bs.se = stat[2,]
	}

	mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.)
	if (anyof(S, 1)) _error(460, "singleton cluster detected")
	if (d==0|d>=.) N = .
	else N = S[,2] :- d

	if (nodots==0) {
		printf("{txt}\nBootstrap replications ({res}%g{txt})\n", reps)
		display("{txt}{hline 4}{c +}{hline 3} 1 " +
		 "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " +
		 "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ")
	}
	bs.rstat = J(reps, cols(bs.stat), .)
	if (hasse) bs.rse = J(reps, cols(bs.se), .)
	for (i=ii=1; i<=reps; i++) {
		ws = mm_sample(N, S, C, 1, 0, 1):*w
		res = mm_callf(setup, X, ws)
		if (missing(res)) {
			if (nodots==0) printf("{err}x{txt}")
		}
		else {
			bs.rstat[ii,] = res[1,]
			if (hasse) bs.rse[ii,] = res[2,]
			ii++
			if (nodots==0) printf(".")
		}
		if (nodots==0) {
			if (mod(i,50)==0) printf(" %5.0f\n",i)
			displayflush()
		}
	}
	if (nodots==0 & mod(i-1,50)) display("")
	if (ii<i) {
		if (ii>1) {
			bs.rstat = bs.rstat[|1,1 \ ii-1,.|]
			if (hasse) bs.rse = bs.rse[|1,1 \ ii-1,.|]
		}
		else {
			bs.rstat = J(0, cols(bs.rstat), .)
			if (hasse) bs.rse = J(0, cols(bs.rse), .)
		}
	}
	bs.reps = reps
	bs.ndrop = reps - rows(bs.rstat)
	return(bs)
}

real matrix mm_bs_report(
 struct mm_bsstats scalar bs, // boostrap estimates
 | string vector what,        // b,mean,bias,v,se,n,basic,p,bc,bca,t
   real scalar level,         //
   real scalar mse,           // !=0 => use mse formula
   struct mm_jkstats jk)      // jackknife estimates
{
	real scalar    theta, mean, bias, v, se, n, basic, p, bc,
	               bca, t, l, i, alpha, tval
	real rowvector bsv, bsse, bsmean, z0, a
	real matrix    res

	if (args()<2) what  = "se"
	if (args()<3) level = st_numscalar("c(level)")
	if (args()<4) mse   = 0

	theta = mean = bias = v = se = n = basic = p = bc = bca = t = l = 0
	for (i=1; i<=length(what); i++) {
		if      (what[i]=="theta" | what[i]=="b") theta = (l = l + 1)
		else if (what[i]=="mean" ) mean  = (l = l + 1)
		else if (what[i]=="bias" ) bias  = (l = l + 1)
		else if (what[i]=="v"    ) v     = (l = l + cols(bs.rstat))
		else if (what[i]=="se"   ) se    = (l = l + 1)
		else if (what[i]=="n" | what[i]=="ci") n = (l = l + 2)
		else if (what[i]=="basic") basic = (l = l + 2)
		else if (what[i]=="p"    ) p     = (l = l + 2)
		else if (what[i]=="bc"   ) bc    = (l = l + 2)
		else if (what[i]=="bca"  ) bca   = (l = l + 2)
		else if (what[i]=="t"    ) t     = (l = l + 2)
		else _error(3498, "'" + what[i] + "' not allowed")
	}

	res = J(l, cols(bs.rstat), .)
// observed parameter vector
	if (theta) res[theta,] = bs.stat
// mean of bs replicates
	if (mean | bias) {
		bsmean = mean(bs.rstat, 1)
		if (mean) res[mean,] = bsmean
	}
// bias
	if (bias) res[bias,] = bsmean - bs.stat
// variance
	if (v) {
		if (mse) bsv = mm_mse(bs.rstat, 1, bs.stat)
		else     bsv = variance(bs.rstat, 1)
		res[|v-cols(res)+1,1 \ v,.|] = bsv
	}
// standard errors
	if (se | n) {
		if (v) bsse = sqrt(diagonal(bsv))'
		else {
			if (mse) bsse = sqrt(mm_colmse(bs.rstat, 1, bs.stat))
			else     bsse = sqrt(mm_colvar(bs.rstat, 1))
		}
		if (se) res[se,] = bsse
	}
// normal ci
	alpha = (100-level)/200
	if (n) {
		tval = invttail(bs.reps-bs.ndrop-1, alpha)
		res[n-1,] = bs.stat - tval*bsse
		res[n,]   = bs.stat + tval*bsse
	}
// basic ci
	if (basic) res[|basic-1,1 \ basic,.|] = 2*bs.stat :-
	 mm_quantile(bs.rstat, 1, ((1-alpha) \ alpha))
// percentile ci
	if (p) res[|p-1,1 \ p,.|] =
	 mm_quantile(bs.rstat, 1, (alpha \ (1-alpha)))
// bias-corrected ci
	if (bc | bca) z0 = editmissing(invnormal(
	 colsum(bs.rstat:<=bs.stat)/rows(bs.rstat) ), 0)
	if (bc) {
		res[|bc-1,1 \ bc,.|] = mm_quantile(bs.rstat, 1,
		 (normal(2*z0 :+ invnormal(alpha)) \
		 normal(2*z0 :- invnormal(alpha))))
	}
// bias-corrected and accelerated ci
	if (bca) {
		a = colsum((mean(jk.rstat, 1):-jk.rstat):^3) :/
		    ( 6 * colsum((mean(jk.rstat, 1):-jk.rstat):^2):^1.5 )
		res[|bca-1,1 \ bca,.|] = mm_quantile(bs.rstat, 1,
		 (normal(z0 :+ (z0 :+ invnormal(alpha)) :/
		   (1:-a:*(z0 :+ invnormal(alpha)))) \
		 normal(z0 :+ (z0 :- invnormal(alpha)) :/
		   (1:-a:*(z0 :- invnormal(alpha))))))
	}
// percentile-t ci
	if (t) {
		res[|t-1,1 \ t,.|] = bs.stat :- bs.se :* mm_quantile(
		 editmissing((bs.rstat:-bs.stat):/bs.rse, 0), 1,
		 ((1-alpha) \ alpha) )
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_rbinomial}{bf:mm_rbinomial.mata}{asis}
*! version 1.0.1, Ben Jann, 14apr2006
version 9.1
mata:

real matrix mm_rbinomial(real matrix n, real matrix p)
{
	real matrix x
	real scalar r, R, c, C
	transmorphic scalar rn, cn, rp, cp

	R = max((rows(n),rows(p)))
	C = max((cols(n),cols(p)))
	rn = (rows(n)==1 ? &1 : (rows(n)<R ? _error(3200) : &r))
	cn = (cols(n)==1 ? &1 : (cols(n)<C ? _error(3200) : &c))
	rp = (rows(p)==1 ? &1 : (rows(p)<R ? _error(3200) : &r))
	cp = (cols(p)==1 ? &1 : (cols(p)<C ? _error(3200) : &c))
	x = J(R,C,.)
	for (r=1;r<=R;r++) {
		for (c=1;c<=C;c++) {
			x[r,c] = _mm_rbinomial(n[*rn,*cn], p[*rp,*cp])
		}
	}
	return(x)
}

real scalar _mm_rbinomial(real scalar n, real scalar p)
{
	if (p<0|p>1) return(.)
	if (n<=0|(n-trunc(n))!=0) return(.)
//rejection technique
	if (n<50|p>.03) {
		return(colsum(uniform(n,1):<p))
	}
//geometric distribution technique
	real scalar sum, x
	sum = 0
	for (x=1;x<=n;x++) {
		sum = sum + ceil((ln(uniform(1,1))/ln(1-p))-1)
		if (sum > n-x) break
	}
	return(x-1)
}

end


*! {smcl}
*! {marker mm_cebinomial}{bf:mm_cebinomial.mata}{asis}
*! version 1.0.1, Ben Jann, 14apr2006
version 9.1
mata:

real matrix mm_cebinomial(real matrix n, real matrix k, real matrix p)
{
	real matrix x
	real scalar r, R, c, C
	transmorphic scalar rn, cn, rk, ck, rp, cp

	R = max((rows(n),rows(k),rows(p)))
	C = max((cols(n),cols(k),cols(p)))
	rn = (rows(n)==1 ? &1 : (rows(n)<R ? _error(3200) : &r))
	cn = (cols(n)==1 ? &1 : (cols(n)<C ? _error(3200) : &c))
	rk = (rows(k)==1 ? &1 : (rows(k)<R ? _error(3200) : &r))
	ck = (cols(k)==1 ? &1 : (cols(k)<C ? _error(3200) : &c))
	rp = (rows(p)==1 ? &1 : (rows(p)<R ? _error(3200) : &r))
	cp = (cols(p)==1 ? &1 : (cols(p)<C ? _error(3200) : &c))
	x = J(R,C,.)
	for (r=1;r<=R;r++) {
		for (c=1;c<=C;c++) {
			x[r,c] = _mm_cebinomial(n[*rn,*cn], k[*rk,*ck], p[*rp,*cp])
		}
	}
	return(x)
}

real scalar _mm_cebinomial(real scalar n, real scalar k, real scalar p)
{
	real scalar e, e0, i

	if (p<0|p>1) return(.)                 //success prob. out of range
	if (n<=0|(n-trunc(n))!=0) return(.)    //n out of range or non-integer
	if (k<0|k>n|(k-trunc(k))!=0) return(.) //k out of range or non-integer
	if (k==0) return(n*p)
	e = 0
	e0 = 0
	for (i=1;i<=n-k;i++) {
		e = e + Binomial(n,k+i,p)
		if (e==e0) break
		e0 = e
	}
	return(e / Binomial(n,k,p) + k)
}
end


*! {smcl}
*! {marker mm_benford}{bf:mm_benford.mata}{asis}
*! version 1.0.0, Ben Jann, 19jun2007
version 9.2
mata:

real vector mm_benford(
 real vector digit0,
 | real scalar pos0,
   real scalar base0)
{
    real scalar    pos, base
    real vector    digit

// check arguments and set defaults
    if (args()<3 | base0>=.) base = 10
    else {
        base = trunc(base0)
        if (base <0 | base >10) _error(3300, "base must be in [2,10]")
    }
    if (args()<2 | pos0>=.) pos = 1
    else {
        pos = trunc(pos0)
        if (pos<1) _error(3300, "pos must be 1 or larger")
    }
    digit = trunc(digit0)
    if (any(digit:<0) | any(digit:>=base))
      _error(3300, "digit must be in [0,base-1] (base is "+strofreal(base)+")")

// case 1: pos==1 (simple formula)
    if (pos==1) return(ln(1 :+ 1:/digit)/ln(base))

// case 2: pos > 1
    if (cols(digit)==1) return( rowsum(ln(1 :+ 1:/(
     J(rows(digit),1,base)*(base^(pos-2)..base^(pos-1)-1) :+ digit))) / ln(base))
    else                return( colsum(ln(1 :+ 1:/(
     (base^(pos-2)::base^(pos-1)-1)*J(1,cols(digit),base) :+ digit))) / ln(base))
}

end


*! {smcl}
*! {marker mm_cauchy}{bf:mm_cauchy.mata}{asis}
version 9.2
// CFBaum 3aug2007
// minor changes by Ben Jann, 09aug2007
// from Wikipedia
// http://en.wikipedia.org/wiki/Cauchy_distribution
mata:
// density function of the Cauchy-Lorentz distribution
real matrix mm_cauchyden(real matrix x, real scalar x0, real scalar gamma)
{
    real matrix cauchyden

    if (gamma <= 0) {
        return(x*.)
    }
    cauchyden = 1 :/ ( pi() * gamma :* ( 1 :+ ( (x :- x0) :/ gamma ) :^2) )
    return(cauchyden)
}

// cumulative distribution function of the Cauchy-Lorentz distribution
real matrix mm_cauchy(real matrix x, real scalar x0, real scalar gamma)
{
    real matrix cauchy

    if (gamma <= 0) {
        return(x*.)
    }
    cauchy = 1 / pi() :* atan( (x :- x0) :/ gamma ) :+ 0.5
    return(cauchy)
}

// tail probability of the Cauchy-Lorentz distribution
// note: cauchy(x,0,1) is the standard Cauchy distribution
//       with cauchytail(x,0,1) = ttail(1,x)
real matrix mm_cauchytail(real matrix x, real scalar x0, real scalar gamma)
{
    return(1 :- mm_cauchy(x, x0, gamma))
}

// inverse cumulative distribution function of the Cauchy-Lorentz distribution
real matrix mm_invcauchy(real matrix p, real scalar x0, real scalar gamma)
{
    real matrix invcauchy

    if ( gamma <= 0 | min(p) < 0 | max(p) > 1 ) {
        return(p*.)
    }
    invcauchy = x0 :+ gamma :* tan( pi() :* ( p :- 0.5 ) ) :+ ln(p:>=0:&p:<=1)
    return(invcauchy)
}
end


*! {smcl}
*! {marker mm_subset}{bf:mm_subset.mata}{asis}
*! version 1.0.1, Ben Jann, 02may2007
version 9.2
mata:

/*-SUBSETS---------------------------------------------------------------*/
// Generate all combinations (subsets) one-by-one (lexicographic)
// based on Algorithm 5.8 from Reingold et al. (1977)
// Usage:
//
//        info = mm_subsetsetup(n,k)
//        while ((c = mm_subset(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// The total number of combinations can be computed using
// Mata's -comb()- function, i.e.
//
//        Ctot = comb(n, k)

struct mm_subsetinfo {
	real scalar    n, k, i, j, counter, algorithm
	real colvector x, y
}

struct mm_subsetinfo scalar mm_subsetsetup(
 real scalar n,
 | real scalar k)
{
	struct mm_subsetinfo scalar s
// setup
	s.n = trunc(n)
	s.k = trunc(k)
	if (s.n>=.) _error(3351)
	if (s.k>=.) s.k = s.n
	else if (s.k>s.n) _error(3300,"k may not be larger than n")
	if (s.n<1|s.k<1) s.i = 0  // => return empty set
	else {
		s.x = -1 \ (1::s.k)
		s.i = 1
	}
	s.i = s.i + 1 // offset i by 1
	s.k = s.k + 1 // offset k by 1
	s.counter = 0
	return(s)
}

real colvector mm_subset(struct mm_subsetinfo scalar s)
{
	real scalar    j
	real colvector subset

	if (s.i==1) return(J(0,1,.)) /*done*/
// get current subset (to be returned at end)
	subset = s.x[|2 \ s.k|]
// generate next subset
	s.i = s.k
	while (s.x[s.i] == s.n-s.k+s.i) {
		s.i = s.i - 1
	}
	s.x[s.i] = s.x[s.i] + 1
	for (j=s.i+1;j<=s.k;j++) s.x[j] = s.x[j-1] + 1
/* the following would improve speed but does not work for some reason
	if (s.i<s.k) s.x[|s.i+1 \ s.k|] = s.x[|s.i \ s.k-1|] :+ 1
*/
	s.counter = s.counter + 1
	return(subset)
}

// wrapper for mm_subsetinfo()/mm_subset()
real matrix mm_subsets(real scalar n, | real scalar k)
{
	real scalar    i
	real colvector set
	real matrix    res
	struct mm_subsetinfo scalar s

	s = mm_subsetsetup(n, (k<. ? k : n))
	if ((set = mm_subset(s)) == J(0,1,.)) return(set)
	res = J(rows(set), comb(n, rows(set)), .)
	res[,i=1] = set
	while ((set = mm_subset(s)) != J(0,1,.)) res[,++i] = set
	return(res)
}


/*-COMPOSITIONS----------------------------------------------------------*/
// Generate all k-way compositions of n, one-by-one
// Default algorithm: direct (less elegant but faster) (anti-lexicographic)
// Alternate algorithm: indirect via subsets (lexicographic)
//
// Usage:
//
//        info = mm_compositionsetup(n,k)
//        while ((c = mm_composition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// or:
//
//        info = mm_compositionsetup(n,k,1)
//        while ((c = mm_composition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// mm_ncompositions() returns the total number of compositions

real scalar mm_ncompositions(
 real scalar n,
 | real scalar k0)
{
	real scalar k

	k = trunc(k0<. ? k0 :  n)
	return(comb(trunc(n)+k-1,k-1))
}

struct mm_subsetinfo scalar mm_compositionsetup(
 real scalar n,
 | real scalar k,
   real scalar alt) // !=0 => alternate algorithm
{
	struct mm_subsetinfo scalar s

	if (args()<3) alt = 0
	if (alt) {
		s = _mm_composition2setup(n,k)
		s.algorithm = 2
		return(s)
	}
	s = _mm_compositionsetup(n,k)
	s.algorithm = 1
	return(s)
}

real colvector mm_composition(struct mm_subsetinfo scalar s)
{
	if (s.algorithm==1) return(_mm_composition(s))
	if (s.algorithm==2) return(_mm_composition2(s))
	_error(3498, "unknown algorithm")
}

// default algorithm

struct mm_subsetinfo scalar _mm_compositionsetup(
 real scalar n,
 | real scalar k)
{
	struct mm_subsetinfo scalar s

	s.n = trunc(n)
	s.k = (k<. ? trunc(k) :  s.n)
	if (s.n>=.) _error(3351)
	if (s.k<1) {
		s.i = 0
		s.x = J(0,1,.)
	}
	else if (s.n<1) {
		s.i = 0
		s.x = J(s.k,1,0)
	}
	else if (s.k<2) {
		s.i = 0
		s.x = J(s.k,1,s.n)
	}
	else {
		s.x = s.n \ J(s.k-1,1,0)
		s.i = 1
	}
	s.j = 2
	s.counter = 0
	return(s)
}

real colvector _mm_composition(struct mm_subsetinfo scalar s)
{
	real colvector subset

	subset = s.x
	if (s.i<1) {
		if (subset!=J(0,1,.)) s.counter = s.counter + 1
		s.x = J(0,1,.)
		return(subset) /*done*/
	}
	while (s.x[s.i]==0 & s.i>1) { // move back
		s.x[s.i] = s.x[s.j]
		s.x[s.j] = 0
		s.j = s.i--
	}
	if (s.x[s.i]>0) {
		s.x[s.i] = s.x[s.i]-1
		s.x[s.j] = s.x[s.j]+1
		if (s.j<s.k) s.i = s.j++
	}
	else {
		s.i = 0
		s.x = J(0,1,.)
	}
	s.counter = s.counter + 1
	return(subset)
}

// alternate algorithm

struct mm_subsetinfo scalar _mm_composition2setup(
 real scalar n,
 | real scalar k0)
{
	real scalar k

	k = trunc(k0<. ? k0 :  n)
	return(mm_subsetsetup(trunc(n)+k-1,k-1))
}

real colvector _mm_composition2(struct mm_subsetinfo scalar s)
{
	real colvector subset

	if (s.i==1) {
		if (s.k==1) {  // length of composition is 1
			s.k = 0
			s.counter = s.counter + 1
			return(s.n)
		}
		return(J(0,1,.)) /*done*/
	}
	subset = mm_subset(s)
	return( (subset \ s.n+1) - (0 \ subset) :- 1 )
}

// wrapper for mm_compositioninfo()/mm_composition()
real matrix mm_compositions(real scalar n, | real scalar k, real scalar alt)
{
	real scalar    i
	real colvector set
	real matrix    res
	struct mm_subsetinfo scalar s

	if (args()<3) alt = 0
	s = mm_compositionsetup(n, (k<. ? k : n), alt)
	if ((set = mm_composition(s)) == J(0,1,.)) return(set)
	res = J(rows(set), mm_ncompositions(n, rows(set)), .)
	res[,i=1] = set
	while ((set = mm_composition(s)) != J(0,1,.)) res[,++i] = set
	return(res)
}

/*-RANDOM SUBSETS--------------------------------------------------------*/
// Return random combination (subsets)
// Usage:
//
//        for (i=1;i<=reps;i++) {
//             c = mm_rsubset(n, k)
//             ... c ...
//        }
//
// mm_rsubset() returns jumbled sets; apply sort() if you need
// ordered sets, e.g.
//
//        c = sort(mm_rsubset(n, k), 1)

real colvector mm_rsubset(
 real scalar n,
 | real scalar k)
{
	if (k<. & k>n) _error(3300, "k may not be larger than n")
	if (n<1|k<1) return(J(0,1,.))
	return(mm_unorder2(n)[|1 \ (k<. ? k : n)|])
}

/*-RANDOM COMPOSITIONS---------------------------------------------------*/
// Return random combination (subsets) and random composition
// Usage:
//
//        for (i=1;i<=reps;i++) {
//             c = mm_rcomposition(n, k)
//             ... c ...
//        }

real colvector mm_rcomposition(
 real scalar n0,
 | real scalar k0)
{
	real scalar    n, k
	real colvector subset

	n = trunc(n0)
	k = (k0<. ? trunc(k0) : n)
	if (n>=.) _error(3351)
	if (k<1) return(J(0,1,.))
	if (n<1) return(J(k,1,0))
	if (k<2) return(J(k,1,n))
	subset = sort(mm_rsubset(n+k-1,k-1), 1)
	return( (subset \ n+k) - (0 \ subset) :- 1 )
}

/*-PARTITIONS------------------------------------------------------------*/
// Generate all k-way partitions of n, one-by-one
// Default algorithm: based on ZS1 algorithm in Zoghbi&Stojmenovic (1998)
//                    (anti-lexicographic)
// Alternate algorithm: based on Algorithm 5.12 in Reingold et al. (1977)
//                    (dictionary order)
// Usage:
//
//        info = mm_partitionsetup(n,k)
//        while ((c = mm_partition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// or:
//
//        info = mm_partitionsetup(n,k,1)
//        while ((c = mm_partition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// mm_npartitions() returns the total number of compositions, based
// on algorithms from
// http://home.att.net/~numericana/answer/numbers.htm#partitions

// number of partitions of n with k or fewer addends (slow)
real scalar mm_npartitions(
 real scalar n,
 | real scalar k)
{
	real scalar      u, i, j, m
	real colvector   p, a

	m = min((n,k))
	if (m==n) return(_mm_npartitions(n))
	p = J(n+1,1,1)
	for (u=2;u<=m;u++) {
		a = p
		p = J(n+1,1,0)
		for (i=0;i<=n;i=i+u) {
			for (j=i+1;j<=n+1;j++) {
				p[j] = p[j] + a[j-i]
			}
		}
	}
	return(p[n+1])
}
// total number of partitions of n (fast)
real scalar _mm_npartitions(
 real scalar n)
{
	real scalar      i, j, s, k
	real colvector   p

	p = J(n+1,1,1)
	for (i=1;i<=n;i++) {
		j = 1 ; k = 1 ; s = 0
		while (j>0) {
			j = i - (3*k*k + k)/2
			if (j>=0) s = s - (-1)^k * p[j+1]
			j = i - (3*k*k - k)/2
			if (j>=0) s = s - (-1)^k * p[j+1]
			k = k + 1
		}
		p[i+1] = s
	}
	return(p[n+1])
}

struct mm_subsetinfo scalar mm_partitionsetup(
 real scalar n,
 | real scalar k,
   real scalar pad,  // zero-padded output
   real scalar alt)  // !=0 => alternate algorithm
{
	struct mm_subsetinfo scalar s

	if (args()<3) pad = 0
	if (args()<4) alt = 0
	if (alt) {
		s = _mm_partition2setup(n,k)
		s.algorithm = 3 + (pad!=0)
		return(s)
	}
	s = _mm_partitionsetup(n,k)
	s.algorithm = 1 + (pad!=0)
	return(s)
}

real colvector mm_partition(struct mm_subsetinfo scalar s)
{
	if (s.algorithm==1) return(_mm_partition(s,0))
	if (s.algorithm==2) return(_mm_partition(s,1))
	if (s.algorithm==3) return(_mm_partition2(s,0))
	if (s.algorithm==4) return(_mm_partition2(s,1))
	_error(3498, "unknown algorithm")
}

// default algorithm

struct mm_subsetinfo scalar _mm_partitionsetup(
 real scalar n,
 | real scalar k)
{
	struct mm_subsetinfo scalar s

	s.n = trunc(n)
	s.k = trunc(k)
	if (s.n>=.) _error(3351)
	if (s.k>=.) s.k = s.n
//	else if (s.k>s.n) _error(3300,"k may not be larger than n")
	if (s.n<1 | s.k<1 ) s.i = -1
	else {
		s.x = s.n \ J(s.k-1, 1, 1)
		s.i = 1
	}
	s.j = 1
	s.counter = 0
	return(s)
}

real colvector _mm_partition(        // original ZS1 Algorithm does not
 struct mm_subsetinfo scalar s,      // support k-way partitions; changes
 real scalar pad)                    // are indicated
{
	real scalar    r, t, l
	real colvector subset

	if (s.i==-1) return(J(0,1,.)) /*done*/
	subset = s.x[|1 \ s.j|]
	if (pad) subset = subset \ J(s.k-rows(subset),1,0)
	l = (s.k<s.n ? s.k : s.n)     // added
	if (s.j>=l) {                 // added; original stopping rule is (s.x[1]==1)
		if (s.x[1]<=s.x[l]+1) {      // changed
			s.i = -1
			s.counter = s.counter + 1
			return(subset)
		}
		while (s.x[s.i]<=s.x[l]+1) {  // step back loop added
			s.j = s.j + s.x[s.i] - 1
			s.i = s.i - 1
		}
	}
	if (s.x[s.i]==2) {
		s.j      = s.j + 1
		s.x[s.i] = 1
		s.x[s.j] = 1             // added
		s.i      = s.i - 1
	}
	else {
		r = s.x[s.i] - 1
		t = s.j - s.i + 1
		s.x[s.i] = r
		while (t>=r) {
			s.i = s.i + 1
			s.x[s.i] = r
			t = t - r
		}
		if (t==0) {
			s.j = s.i
		}
		else {
			s.j = s.i + 1
			if (t==1) s.x[s.i+1] = 1    // added
			if (t>1) {
				s.i    = s.i + 1
				s.x[s.i] = t
			}
		}
	}
	s.counter = s.counter + 1
	return(subset)
}

// alternate algorithm

struct mm_subsetinfo scalar _mm_partition2setup(
 real scalar n,
 | real scalar k)
{
	real scalar offset
	struct mm_subsetinfo scalar s

	s.n = trunc(n)
	s.k = trunc(k)
	if (s.n>=.) _error(3351)
	if (s.k>=.) s.k = s.n
	offset = 2
	if (s.n<1 | s.k<1 ) s.i = 0 + offset
	else {
		s.i = 1 + offset
		s.x = s.y = J(s.n+offset, 1, .)
		s.x[-1+offset] = 0
		s.y[-1+offset] = 0
		s.x[0+offset]  = 0
		s.y[0+offset]  = s.n + 1
		s.x[1+offset]  = s.n
		s.y[1+offset]  = 1
	}
	s.counter = 0
	return(s)
}

real colvector _mm_partition2(
 struct mm_subsetinfo scalar s,
 real scalar pad)
{
	real scalar    skip, sum
	real colvector subset

	skip = 1
	while (skip) {
		if (s.i<=2) return(J(0,1,.)) /*done*/
		if ((skip = (sum(s.x[|3 \ s.i|])>s.k))==0)
		 subset = _mm_partition2_expand(s.x[|3 \ s.i|], s.y[|3 \ s.i|],
		 (pad ? s.k : .))
		sum = s.x[s.i]*s.y[s.i]
		if (s.x[s.i]==1) {
			s.i = s.i - 1
			sum = sum + s.x[s.i]*s.y[s.i]
		}
		if (s.y[s.i-1]==s.y[s.i]+1) {
			s.i      = s.i - 1
			s.x[s.i] = s.x[s.i] + 1
		}
		else {
			s.y[s.i] = s.y[s.i] + 1
			s.x[s.i] = 1
		}
		if (sum>s.y[s.i]) {
			s.y[s.i+1] = 1
			s.x[s.i+1] = sum - s.y[s.i]
			s.i        = s.i + 1
		}
	}
	s.counter = s.counter + 1
	return(subset)
}

real colvector _mm_partition2_expand(
 real vector m,
 real vector p,
 real scalar k)
{
	real scalar    i, j, l
	real colvector res

	if (k<.) res = J(k,1,0)
	else     res = J(sum(m),1,0)
	l = 0
	for (i=1;i<=rows(p);i++) {
		for (j=1;j<=m[i];j++) {
			res[++l] = p[i]
		}
	}
	return(res)
}

// wrapper for mm_partitioninfo()/mm_partition()
real matrix mm_partitions(real scalar n, | real scalar k, real scalar alt)
{
	real scalar    i
	real colvector set
	real matrix    res
	struct mm_subsetinfo scalar s

	if (args()<3) alt = 0
	s = mm_partitionsetup(n, (k<. ? k : n), 1, alt)
	if ((set = mm_partition(s)) == J(0,1,.)) return(set)
	res = J(rows(set), mm_npartitions(n, rows(set)), .)
	res[,i=1] = set
	while ((set = mm_partition(s)) != J(0,1,.)) res[,++i] = set
	return(res)
}

end


*! {smcl}
*! {marker mm_mgof}{bf:mm_mgof.mata}{asis}
*! version 1.0.7  Ben Jann  12jan2008
*  mm_mgof: goodness-of-fit tests for multinomial data

version 9.2
mata:

real matrix mm_mgof(
    real colvector f,    // observed count
  | real colvector h0,   // expected distribution
    string scalar m0,    // method: "approx" (default), "mc", "rc", or "ee"
    string vector s0,    // test-stats: "x2", "lr", "cr", "mlnp", "ksmirnov"
    real scalar lambda,  // lambda for Cressie-Read statistic
    real scalar arg1,    // nfit for "approx"; reps for "mc"
    real scalar dots,    // display dots with mc/ee
    real scalar arg2)    // will be replaced; used by mgof_ee to return counter
{
    real scalar     i, n, hsum
    string scalar   m, s
    string vector   stats
    real vector     p
    pointer(real colvector) scalar       h
    pointer(real scalar function) vector s1

// parsing and defaults
    if (args()<7) dots = 0
    m = mm_strexpand(m0, ("approx", "mc", "rc", "ee"), "approx", 0)
    stats = ("x2", "lr", "cr")
    if (m!="approx") stats = (stats, "mlnp", "ksmirnov")
    if (length(s0)<1) s = "x2"
    else {
        s = J(length(s0),1,"")
        for (i=1;i<=length(s0);i++) {
            s[i] = mm_strexpand(s0[i], stats, "x2", 1)
        }
    }
    if (m!="approx" & lambda<0) _error(3498, "lambda<0 not allowed with "+m)

// check data / prepare null distribution
    if (missing(f)) _error(3351, "f has missing values")
    if (any(f:<0)) _error(3498, "f has negative values")
    if (any(f:<0)) _error(3498, "f has negative values")
    if (m=="ee" | m=="rc" | (m=="mc" & anyof(s, "mlnp"))) {
        if (f!=trunc(f)) _error(3498, "non-integer f not allowed in this context")
    }
    if (anyof(stats,"cr") & lambda<0 & anyof(f,0))
     _error(3498, "some f are 0: lambda<0 not allowed")
    if (args()<2) h0 = 1
    if (rows(f)!=rows(h0) & rows(h0)!=1) _error(3200)
    if (missing(h0)) _error(3351, "h has missing values")
    if (any(h0:<=0)) _error(3498, "some h are negative or 0")
    n = colsum(f)
    hsum = colsum(h0)
    if (n!=hsum) h = &(h0*n/hsum)
    else h = &h0

// return (0,1) if n. of obs is 0 or n. of categories < 2
    if ( n==0 | rows(f)<2 ) return(J(rows(s),1,0),J(rows(s),1,1))

// gather statistics functions
    s1 = J(rows(s),1,NULL)
    for (i=1;i<=length(s);i++) {
        if      (s[i]=="x2")        s1[i] = &_mm_mgof_x2()
        else if (s[i]=="lr")        s1[i] = &_mm_mgof_lr()
        else if (s[i]=="cr")        s1[i] = &_mm_mgof_cr()
        else if (s[i]=="mlnp")      s1[i] = &_mm_mgof_mlnp()
        else if (s[i]=="ksmirnov")  s1[i] = &_mm_mgof_ksmirnov()
        else _error(3498,s[i]+ " invalid")
    }

// approximate chi2-test
    if (m=="approx") return(_mm_mgof_approx(s1, f, *h, lambda, arg1))

// Monte Carlo test
    if (m=="mc") return(_mm_mgof_mc(s1, f, *h, lambda, arg1, dots))

// exhaustive enumeration / random composition test
    if (s=="mlnp") {
        if (m=="rc") return(_mm_mgof_rc(NULL, f, *h))
        return(_mm_mgof_ee(NULL, f, *h))
    }
    p = mm_which(s:=="mlnp")
    if (rows(p)>0) {
        p = p[1]
        p = p \ (p>1 ? (1::p-1) : J(0,1,.)) \ (p<rows(s) ? (p+1::rows(s)) : J(0,1,.))
        s1 = s1[p[| 2 \ .|]]
        p = invorder(p)
    }
    else p = (2::length(s)+1)
    if (m=="rc") return(_mm_mgof_rc(s1, f, *h, lambda, arg1)[p,])
    return(_mm_mgof_ee(s1, f, *h, lambda, anyof(s,"ksmirnov"), dots, arg2)[p,])
}

// large sample chi2-approximation test for multinomial distributions
real matrix _mm_mgof_approx(
 pointer(real scalar function) vector s, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar nfit)   // n. of fitted parameters (default 0)
{
    real scalar                     n, i, k, nstat
    real colvector                  stat
    pointer(real colvector) scalar  h

// compute hypothetical distribution if missing
    n = colsum(f)
    k = rows(f)
    if (rows(h0)<=1) h = &(n :* J(k,1,1/k))
    else             h = &h0
// compute statistics and p-values
    nstat = length(s)
    stat = J(nstat,1,.)
    for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h,lambda)
    return(stat, chi2tail(rows(f) - (nfit<. ? nfit : 0) - 1,stat))
}

// Monte Carlo exact test for discrete distributions
// (based on sampling from the hypothetical distribution)
real matrix _mm_mgof_mc(
 pointer(real scalar function) vector s, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar reps0,  // replications (default=10000)
   real scalar dots)   // display dots
{
    real scalar                     n, j, i, k, nstat, reps, uni
    real scalar                     doti, ndots
    real colvector                  fj, H, stat, p, hp
    pointer(real colvector) scalar  h

// set defaults and compute hypothetical distribution
    if (args()<5 | reps0>=.) reps = 10000
    else                     reps = reps0
    if (args()<6) dots = 0
    n = round(colsum(f)) // so that, e.g., 10 are sampled if n=9.9999...
    k = rows(f)
    if (uni = (rows(h0)==1)) {
        uni = 1
        h  = &(n :* J(k,1,1/k))
        H  = rangen(1/k,1,k)
        hp = J(k,1,1/k)
    }
    else {
        uni = mm_isconstant(h0)
        H  = mm_colrunsum(h0)
        H  = H / H[rows(H)]
        hp = (h0) / colsum(h0)
        h  = &h0
    }
// compute statistics and p-values
    nstat = length(s)
    stat = p = J(nstat,1,0)
    for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H)
    if (reps<=0) return(stat, J(nstat,1,.))
    if (dots) {
        printf("\n{txt}Percent completed ({res}%g{txt} replications)\n",reps)
        display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hline 5} 100")
        ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 :
                (reps>=5 ? 10 : (reps>=2 ? 25 : 50)))))
        doti = 1
    }
    for (j=1; j<=reps; j++) {
        if (uni) fj = _mm_mgof_mc_usmpl(n, k) // = mm_srswr(n,rows(f),1)
        else     fj = _mm_mgof_mc_smpl(n, H) // = mm_upswr(n, h, 1)  (slower)
        for (i=1; i<=nstat; i++)  p[i] = p[i] +
          (stat[i] <= (*s[i])(fj,*h, lambda, n, hp, H))
        if (dots) {
            if (j*50 >= doti*ndots*reps) {
                printf(ndots*".")
                displayflush()
                doti++
            }
        }
    }
    if (dots) display("")
    return(stat, p/reps)
}
// Monte Carlo subroutine: sample from uniform multinomial
real colvector _mm_mgof_mc_usmpl(
 real scalar n, // sample size
 real scalar k) // number of categories
{
    real scalar     i
    real colvector  u, res

    u = ceil(uniform(n,1)*k)
    res = J(k,1,0)
    if (n>1.5^(k-2)) { // faster for small k relative to n
        for (i=1; i<=k; i++) res[i] = colsum(u:==i)
    }
    else {
        for (i=1;i<=n;i++) res[u[i]] = res[u[i]] + 1
    }
    return(res)
}
// Monte Carlo subroutine: sample from non-uniform multinomial
real colvector _mm_mgof_mc_smpl(
 real scalar n,    // sample size
 real colvector F) // distribution function (cumulative)
{
    real scalar     i, j
    real colvector  u, res

    if (n>2.3^(rows(F)-3)) {  // faster for small k relative to n
        u = uniform(n,1)
        res = J(rows(F),1,0)
        res[1] = sum(u:<=F[1])
        for (j=2;j<=rows(F);j++) {
            res[j] = sum(u:>F[j-1] :& u:<=F[j])
        }
        return(res)
    }
    u = sort(uniform(n,1),1)
    j=1
    res = J(rows(F),1,0)
    for (i=1;i<=n;i++) {
        while (u[i]>F[j]) j++
        res[j] = res[j]+1
    }
    return(res)
}

// Random composition exact test for discrete distributions
// -> undocumented; do not use this
//    the method is highly unreliable if the number of potential
//    compositions is large
real matrix _mm_mgof_rc(
 pointer(real scalar function) vector s0, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar reps0)  // replications (default=10000)
{
    real scalar                           n, i, j, k, nstat, reps, lnscale
    real colvector                        fj, H, stat, statj, p, jj, se, hp
    pointer(real colvector) scalar        h
    pointer(real scalar function) vector  s

// set defaults and compute hypothetical distribution
    if (args()<5 | reps0>=.) reps = 10000
    else                     reps = reps0
    if (s0==NULL) s = &_mm_mgof_mlnp()
    else {
        if (rows(s0)!=1) s = &_mm_mgof_mlnp() \ s0
        else             s = &_mm_mgof_mlnp() , s0
    }
    n = colsum(f)
    k = rows(f)
    if (rows(h0)==1) {
        h = &(n :* J(k,1,1/k))
        H = rangen(1/k,1,k)
        hp = J(k,1,1/k)
    }
    else {
        h = &h0
        H = mm_colrunsum(*h)
        H = H / H[rows(H)]
        hp = (*h) / colsum(*h)
    }
// compute statistics and p-values
    nstat = length(s)
    stat = statj = p = jj = se = J(nstat,1,0)
    for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H)
    if (reps<=0) return(stat, J(nstat,2,.))
    for (j=1; j<=reps; j++) {
        fj = mm_rcomposition(n, k)
        for (i=1; i<=nstat; i++) {
            statj[i] = (*s[i])(fj,*h, lambda, n, hp, H)
            if (stat[i] <= statj[i]) {
                jj[i] = jj[i] + 1
                p[i]  = p[i] + (exp(-statj[1])-p[i])/jj[i]     // mean updating
                se[i] = se[i] + (exp(-statj[1])^2-se[i])/jj[i] // mean updating
            }
        }
    }
    p = p:*jj; se = se:*jj
    // lnscale = log. of inverse sampling probability
    lnscale = lnfactorial(n+k-1) - lnfactorial(k-1) - lnfactorial(n) - ln(reps)
    se = exp( lnscale :+ ln(sqrt(reps*(1/(reps-1) * se - (p/reps):^2))) )
//  se = mm_ncompositions(n, k) * sqrt((1/(reps-1) * se - (p/reps):^2)/reps)
    p  = exp( lnscale :+ ln(p) )
//  p  = mm_ncompositions(n, k) / reps * p
    return(stat, p, se)
}

// exhaustive enumeration exact test (performs the exact
// multinomial g.o.f. test plus, optionally, exact tests
// for specified additional statistics)
// WARNING: use this test for very small samples only since
// computation time is linear to the number of compositions
// =comb(n+k-1,k-1), where n is the sample size and k is
// the number of categories.
// IMPORTANT EXCEPTION: when the null distribution is uniform,
// compution time is determined by the number of partitions
// which is magnitudes smaller than the number of compositions
real matrix _mm_mgof_ee(
 pointer(real scalar function) vector s0, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar force,  // do not use partitions (for ksmirnov)
   real scalar dots,   // display dots
   real scalar counter) // will be replaced
{
    real scalar                           n, i, k, nstat, uni, w
    real scalar                           doti, ndots, reps
    real colvector                        fj, j, H, stat, statj, p, hp
    pointer(real colvector) scalar        h
    pointer(real scalar function) vector  s
    struct mm_subsetinfo scalar           info

// set defaults and compute hypothetical distribution
    if (args()<6) dots = 0
    if (s0==NULL) s = &_mm_mgof_mlnp()
    else {
        if (rows(s0)!=1) s = &_mm_mgof_mlnp() \ s0
        else             s = &_mm_mgof_mlnp() , s0
    }
    n = colsum(f)
    k = rows(f)
    uni = (force==1 ? 0 : 1)
    if (rows(h0)==1) {
        h = &(n :* J(k,1,1/k))
        H = rangen(1/k,1,k)
        hp = J(k,1,1/k)
    }
    else {
        h = &h0
        if (mm_isconstant(*h)==0) uni = 0 // reset to 0 if h not constant
        H = mm_colrunsum(*h)
        H = H / H[rows(H)]
        hp = (*h) / colsum(*h)
    }
// compute statistics and p-values
    nstat = length(s)
    stat = statj = p = j = J(nstat,1,0)
    for (i=1; i<=nstat; i++)  stat[i] = (*s[i])(f,*h, lambda, n, hp, H)
    if (uni) {         // uniform null distribution: work with partitions
        if (dots) {
            reps = mm_npartitions(n,k)
            printf("\n{txt}Percent completed ({res}%g{txt} partitions)\n",reps)
            display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hline 5} 100")
            ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 :
                    (reps>=5 ? 10 : (reps>=2 ? 25 : 50)))))
            doti = 1
        }
        info = _mm_partitionsetup(n, k)
        while ((fj = _mm_partition(info,1)) != J(0,1,.)) {
            w = exp(lnfactorial(info.k) - quadcolsum(lnfactorial(_mm_panels(fj))))
            for (i=1; i<=nstat; i++) {
                statj[i] = (*s[i])(fj,*h, lambda, n, hp, H)
                if (stat[i] <= statj[i]) {
                    j[i] = j[i] + w // statj[1] = -ln(p) => p = exp(-statj[1])
                    p[i] = p[i] + (exp(-statj[1])-p[i])*w/j[i]  // mean updating
                }
            }
            if (dots) {
                if (info.counter*50 >= doti*ndots*reps) {
                    printf(ndots*".")
                    displayflush()
                    doti++
                }
            }
        }
        if (dots) display("")
        counter = info.counter
        return(stat, p:*j)
    }
    if (dots) {
        reps = mm_ncompositions(n,k)
        printf("\n{txt}Percent completed ({res}%g{txt} compositions)\n",reps)
        display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hline 5} 100")
        ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 :
                (reps>=5 ? 10 : (reps>=2 ? 25 : 50)))))
        doti = 1
    }
    info = _mm_compositionsetup(n, k)
    while ((fj = _mm_composition(info)) != J(0,1,.)) {
        for (i=1; i<=nstat; i++) {
            statj[i] = (*s[i])(fj,*h, lambda, n, hp, H)
            if (stat[i] <= statj[i]) {
                j[i] = j[i] + 1  // statj[1] = -ln(p) => p = exp(-statj[1])
                p[i] = p[i] + (exp(-statj[1])-p[i])/j[i]  // mean updating
            }
        }
        if (dots) {
            if (info.counter*50 >= doti*ndots*reps) {
                printf(ndots*".")
                displayflush()
                doti++
            }
        }
    }
    if (dots) display("")
    counter = info.counter
    return(stat, p:*j)
}

// Pearson's chi2 statistic subroutine
real scalar _mm_mgof_x2(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      transmorphic matrix opt2, // not used
      transmorphic matrix opt3, // not used
      transmorphic matrix opt4) // not used
{
    return(quadcolsum((f-h):^2:/h))
}

// log likelihood ratio statistic subroutine
real scalar _mm_mgof_lr(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      transmorphic matrix opt2, // not used
      transmorphic matrix opt3, // not used
      transmorphic matrix opt4) // not used
{
    return(2*quadcolsum(f:*(ln(f)-ln(h))))
}

// Cressie-Read statistic
real scalar _mm_mgof_cr(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | real scalar lambda,       // lambda parameter; default is 2/3
      transmorphic matrix opt2, // not used
      transmorphic matrix opt3, // not used
      transmorphic matrix opt4) // not used
{
    real scalar L
    L = (lambda<. ? lambda : 2/3)
    if (abs(L)<1e-6)        return( 2 * quadcolsum(f:*(ln(f)-ln(h))) )
    else if (abs(L+1)<1e-6) return( 2 * quadcolsum(h:*(ln(h)-ln(f))) )
    else if (L>0) return( 2/(L*(L+1)) * quadcolsum(f:*((f:/h):^L :- 1)) )
                  return( 2/(L*(L+1)) * quadcolsum(f:*((h:/f):^-L :- 1)) )
}

// -ln(p) statistic subroutine (multinomial test)
real scalar _mm_mgof_mlnp(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      real scalar n,            // n. of obs [ =colsum(f) ]
      real colvector p,         // hypothical propotions [ = h/colsum(h) ]
      transmorphic matrix opt4) // not used
{
    if (args()>=5) return(-(lnfactorial(n) -
     quadcolsum(lnfactorial(f)) + quadcolsum(f:*ln(p))))
    return(-(lnfactorial(quadcolsum(f)) -
     quadcolsum(lnfactorial(f)) + quadcolsum(f:*ln(h/quadcolsum(h)))))
}

// kolmogorov-smirnov statistic subroutine
real scalar _mm_mgof_ksmirnov(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      real scalar n,            // n. of obs [ =colsum(f) ]
      transmorphic matrix opt3, // not used
      real colvector H)         // hypothetical cumulative [ = runsum(h)/sum(h) ]
{
    if (args()>=6) return(max(abs(H-mm_colrunsum(f)/n)))
    return(max(abs(mm_colrunsum(h)/colsum(h)-mm_colrunsum(f)/colsum(f))))
}


end


*! {smcl}
*! {marker mm_root}{bf:mm_root.mata}{asis}
*! version 1.0.0, Ben Jann, 07jul2006
*! translation of zeroin.c from http://www.netlib.org/c/brent.shar

* Quoted from http://www.netlib.org/c/brent.shar:
* X * Algorithm
* X *    G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
* X *    computations. M., Mir, 1980, p.180 of the Russian edition
* X *
* X *    The function makes use of the bissection procedure combined with
* X *    the linear or quadric inverse interpolation.
* X *    At every step program operates on three abscissae - a, b, and c.
* X *    b - the last and the best approximation to the root
* X *    a - the last but one approximation
* X *    c - the last but one or even earlier approximation than a that
* X *        1) |f(b)| <= |f(c)|
* X *        2) f(b) and f(c) have opposite signs, i.e. b and c confine
* X *           the root
* X *    At every step Zeroin selects one of the two new approximations, the
* X *    former being obtained by the bissection procedure and the latter
* X *    resulting in the interpolation (if a,b, and c are all different
* X *    the quadric interpolation is utilized, otherwise the linear one).
* X *    If the latter (i.e. obtained by the interpolation) point is
* X *    reasonable (i.e. lies within the current interval [b,c] not being
* X *    too close to the boundaries) it is accepted. The bissection result
* X *    is used in the other case. Therefore, the range of uncertainty is
* X *    ensured to be reduced at least by the factor 1.6

* Main changes compared to the original:
* 1. mm_root() returns an error code and stored the solution in x.
* 2. The maximum # of iteration may be set.
* 3. The tolerance argument is optional.
* 4. Additional arguments may be passed on to f.
* 5. Special action is taken if f returns missing, convergence is not
*    reached, or f(ax) and f(bx) do not have opposite signs.

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real scalar mm_root(
 transmorphic x,      // bj: will be replaced by solution
 pointer(real scalar function) scalar f,
                      // Address of the function whose zero will be sought for
 real scalar ax,      // Root will be sought for within a range [ax,bx]
 real scalar bx,      //
 | real scalar tol,   // Acceptable tolerance for the root value (default 0)
   real scalar maxit, // bj: maximum # of iterations (default: 1000)
   `opts')            // bj: additional args to pass on to f
{
    transmorphic  fs            // setup for f
    real scalar   a, b, c       // Abscissae, descr. see above
    real scalar   fa, fb, fc    // f(a), f(b), f(c)
    real scalar   prev_step     // Distance from the last but one
    real scalar   tol_act       // Actual tolerance
    real scalar   p             // Interpolation step is calcu-
    real scalar   q             // lated in the form p/q; divi-
                                // sion operations is delayed
                                // until the last moment
    real scalar   new_step      // Step at this iteration
    real scalar   t1, cb, t2
    real scalar   itr

    if (args()<5) tol = 0       // bj: set tolerance
    if (args()<6) maxit = 1000  // bj: maximum # of iterations

    fs = mm_callf_setup(f, args()-6, `opts') // bj: prepare function call

    x = .                       // bj: initialize output

    a = ax;  b = bx;  fa = mm_callf(fs, a);  fb = mm_callf(fs, b)
    c = a;  fc = fa

    if ( fa==. ) return(0)      // bj: abort if fa missing

    if ( (fa > 0 & fb > 0) |    // bj: f(ax) and f(bx) do not
         (fa < 0 & fb < 0) ) {  //     have opposite signs
        if ( abs(fa) < abs(fb) ) {
            x = a; return(2)    // bj: fa closer to zero than fb
        }
        x = b; return(3)        // bj: fb closer to zero than fa
    }

    for (itr=1; itr<=maxit; itr++) {
        if ( fb==. ) return(0)            // bj: abort if fb missing

        prev_step = b-a

        if( abs(fc) < abs(fb) ) {         // Swap data for b to be the
            a = b;  b = c;  c = a;        // best approximation
            fa = fb;  fb = fc;  fc = fa
        }

        tol_act = 2*epsilon(b) + tol/2
        new_step = (c-b)/2

        if( abs(new_step) <= tol_act | fb == 0 ) {
             x = b                        // Acceptable approx. is found
             return(0)
        }

        // Decide if the interpolation can be tried
        if( abs(prev_step) >= tol_act     // If prev_step was large enough
             & abs(fa) > abs(fb) ) {      // and was in true direction,
                                          // Interpolation may be tried
            cb = c-b
            if( a==c ) {                  // If we have only two distinct
                t1 = fb/fa                // points linear interpolation
                p = cb*t1                 // can only be applied
                q = 1.0 - t1
            }
            else {                        // Quadric inverse interpolation
                q = fa/fc;  t1 = fb/fc;  t2 = fb/fa
                p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) )
                q = (q-1.0) * (t1-1.0) * (t2-1.0)
            }
            if( p>0 )                     // p was calculated with the op-
              q = -q                      // posite sign; make p positive
            else                          // and assign possible minus to
              p = -p                      // q
            if( p < (0.75*cb*q-abs(tol_act*q)/2) // If b+p/q falls in [b,c]
               & p < abs(prev_step*q/2) ) // and isn't too large
             new_step = p/q               // it is accepted
                                          // If p/q is too large then the
                                          // bisection procedure can
                                          // reduce [b,c] range to more
                                          // extent
        }

        if( abs(new_step) < tol_act ) {   // Adjust the step to be not less
            if( new_step > 0 )            // than tolerance
             new_step = tol_act
            else
             new_step = -tol_act
        }

        a = b;  fa = fb                   // Save the previous approx.
        b = b + new_step;  fb = mm_callf(fs, b) // Do step to a new approxim.
        if( (fb > 0 & fc > 0) | (fb < 0 & fc < 0) )  {
                      // Adjust c for it to have a sign opposite to that of b
            c = a;  fc = fa
        }
    }
    x = b
    return(1)                             // bj: convergence not reached
}

end


*! {smcl}
*! {marker mm_nrroot}{bf:mm_nrroot.mata}{asis}
*! version 1.0.0, Ben Jann, 07jul2006

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real scalar mm_nrroot(
 real scalar x,              // initial guess; will be replaced by solution
 pointer(function) scalar f, // the function; must return f() and f'()
 | real scalar tol,          // acceptable tolerance (default 0)
   real scalar maxit,        // maximum # of iterations (default: 1000)
   `opts')                   // additional args to pass on to f
{
    transmorphic  fs         // setup for f
    real scalar   itr
    real scalar   fx         // fx = ( f(x), f'(x) )
    real scalar   dx         // dx = fx[1] / fx[2]
    real scalar   tol_act    // actual tolerance

    if (args()<3) tol = 0      // default tolerance
    if (args()<4) maxit = 1000 // default maximum # of iterations

    fs = mm_callf_setup(f, args()-4, `opts') // prepare function call
    for (itr=1; itr<=maxit; itr++) {
        tol_act = 2e+4*epsilon(x) + tol/2
        fx = mm_callf(fs, x)
        dx = fx[1] / fx[2]
        x = x - dx
        if ( abs(dx) <= tol_act | missing(x) ) return(0)
    }
    return(1)                // convergence not reached
}

end


*! {smcl}
*! {marker mm_minim}{bf:mm_minim.mata}{asis}
*! version 1.0.0  Ben Jann  05aug2020
*! translation of function Brent_fmin() from file optimize.c included in the
*! source of R 4.0.2

/*
Quoted from optimize.c:
    an approximation  x  to the point where  f  attains a minimum  on
the interval  (ax,bx)  is determined.

INPUT..

ax    left endpoint of initial interval
bx    right endpoint of initial interval
f     function which evaluates  f(x, info)  for any  x
      in the interval  (ax,bx)
tol   desired length of the interval of uncertainty of the final
      result ( >= 0.)

OUTPUT..

fmin  abcissa approximating the point where  f  attains a minimum

    The method used is a combination of  golden  section  search  and
successive parabolic interpolation.  convergence is never much slower
than  that  for  a  Fibonacci search.  If  f  has a continuous second
derivative which is positive at the minimum (which is not  at  ax  or
bx),  then  convergence  is  superlinear, and usually of the order of
about  1.324....
    The function  f  is never evaluated at two points closer together
than  eps*abs(fmin)+(tol/3), where eps is  approximately  the  square
root  of  the  relative  machine  precision.   if   f   is a unimodal
function and the computed values of   f   are  always  unimodal  when
separated  by  at least  eps*abs(x)+(tol/3), then  fmin  approximates
the abcissa of the global minimum of  f  on the interval  ax,bx  with
an error less than  3*eps*abs(fmin)+tol.  if   f   is  not  unimodal,
then fmin may approximate a local, but perhaps non-global, minimum to
the same accuracy.
    This function subprogram is a slightly modified  version  of  the
Algol  60 procedure  localmin  given in Richard Brent, Algorithms for
Minimization without Derivatives, Prentice-Hall, Inc. (1973).
*/

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real scalar mm_minim(
    pointer(real scalar function) scalar f, // address of function to be minimized
    real scalar ax,   // minimum will be sought for within [ax,bx]
    real scalar bx,
  | real scalar tol0, // acceptable tolerance (length of uncertainty interval)
    `opts')           // additional args to pass on to f
{
    transmorphic  fs            // setup for f
    real scalar   a, b, c, d, e, p, q, r, u, v, w, x
    real scalar   t2, fu, fv, fw, fx, xm, eps, tol1, tol3, tol
    
    if (tol0<=0) _error(3300)
    tol = (tol0>=. ? epsilon(1)^0.25 : tol0) // default tolerance
    fs = mm_callf_setup(f, args()-4, `opts') // prepare function call
    c = (3 - sqrt(5)) * .5                   // squared inverse of the golden ratio
    eps = epsilon(1)
    tol1 = eps + 1 // the smallest 1.000... > 1
    eps = sqrt(eps)
    a = ax; b = bx
    v = w = x = a + c * (b - a)
    d = e = 0
    fv = fw = fx = mm_callf(fs, x)
    tol3 = tol / 3
    while (1) {
        xm = (a + b) * .5
        tol1 = eps * abs(x) + tol3
        t2 = tol1 * 2
        /* check stopping criterion */
        if (abs(x - xm) <= t2 - (b - a) * .5) break
        p = q = r = 0
        /* fit parabola */
        if (abs(e) > tol1) {
            r = (x - w) * (fx - fv)
            q = (x - v) * (fx - fw)
            p = (x - v) * q - (x - w) * r
            q = (q - r) * 2
            if (q > 0) p = -p
            else       q = -q
            r = e
            e = d
        }
        /* a golden-section step */
        if (abs(p) >= abs(q * .5 * r) | p <= q * (a - x) | p >= q * (b - x)) {
            if (x < xm) e = b - x
            else        e = a - x
            d = c * e
        }
        /* a parabolic-interpolation step */
        else {
            d = p / q
            u = x + d
            /* f must not be evaluated too close to ax or bx */
            if (u - a < t2 | b - u < t2) {
                d = tol1
                if (x >= xm) d = -d
            }
        }
        /* f must not be evaluated too close to x */
        if (abs(d) >= tol1) u = x + d
        else if (d > 0)     u = x + tol1
        else                u = x - tol1
        fu = mm_callf(fs, u)
        /*  update  a, b, v, w, and x */
        if (fu <= fx) {
            if (u < x) b = x
            else       a = x
            v = w;    w = x;   x = u
            fv = fw; fw = fx; fx = fu
        }
        else {
            if (u < x) a = u
            else       b = u
            if (fu <= fw | w == x) {
                v = w; fv = fw
                w = u; fw = fu
            } 
            else if (fu <= fv | v == x | v == w) {
                v = u; fv = fu
            }
        }
    }
    return(x)
}

end


*! {smcl}
*! {marker mm_finvert}{bf:mm_finvert.mata}{asis}
*! version 1.0.1, Ben Jann, 16may2014

version 9.2

mata:

real vector mm_finvert(
    y,
    f,     //     nr:      brent:
    o1,    //     fd       lo
    | o2,  //     x0       up
      tol,
      maxit,
      opt)
{
    if (args()<5) tol = 0
    if (args()<6) maxit = 1000
    if (ispointer(o1)) {
        if (args()<4) o2 = y   // set initial guess to y
        if (args()<7)
            return(_mm_finvert_nr(y, f, o1, o2, tol, maxit))
        else
            return(_mm_finvert_nr(y, f, o1, o2, tol, maxit, opt))
    }
    if (args()==3)
        _error(3001, "expected 4 to 6 arguments but received 3")
    if (args()<7)
        return(_mm_finvert_brent(y, f, o1, o2, tol, maxit))
    else 
        return(_mm_finvert_brent(y, f, o1, o2, tol, maxit, opt))
}

real vector _mm_finvert_nr(
    real vector y,
    pointer(real scalar function) scalar f,
    pointer(real scalar function) scalar df,
    real vector x,
    real scalar tol,
    real scalar maxit,
    | opt)
{
    real scalar     i, I, resi, rc
    real vector     res
    pointer scalar  ix

    I = length(y)
    ix = (length(x)==1 ? &1 : (length(x)<I ? _error(3200) : &i))

    res = J(rows(y), cols(y), .)
    for (i=1; i<=I; i++) {
        if (args()<7)
            rc = mm_nrroot(resi=x[*ix], &_mm_finvert_nr_fdf(),
                tol, maxit, y[i], f, df)
        else
            rc = mm_nrroot(resi=x[*ix], &_mm_finvert_nr_fdf(),
                tol, maxit, y[i], f, df, opt)
        if (rc)
         _error(3360, "failure to converge for element " + strofreal(i))
        res[i] = resi
    }
    return(res)
}

real rowvector _mm_finvert_nr_fdf(
    real scalar x,
    real scalar y,
    pointer(real scalar function) scalar f,
    pointer(real scalar function) scalar df,
    | opt)
{
    if (args()<5)
        return( (*f)(x) - y, (*df)(x) )
    else
        return( (*f)(x, opt) - y, (*df)(x, opt) )
}

real vector _mm_finvert_brent(
    real vector y,
    pointer(real scalar function) scalar f,
    real vector lo,
    real vector up,
    real scalar tol,
    real scalar maxit,
    | opt)
{
    real scalar     i, I, resi, rc
    real vector     res
    pointer scalar  il, iu

    I = length(y)
    il = (length(lo)==1 ? &1 : (length(lo)<I ? _error(3200) : &i))
    iu = (length(up)==1 ? &1 : (length(up)<I ? _error(3200) : &i))

    res = J(rows(y), cols(y), .)
    for (i=1; i<=I; i++) {
        if (args()<7)
            rc = mm_root(resi=., &_mm_finvert_brent_f(),
                lo[*il], up[*iu], tol, maxit, y[i], f)
        else
            rc = mm_root(resi=., &_mm_finvert_brent_f(),
                lo[*il], up[*iu], tol, maxit, y[i], f, opt)
        if (rc==1)
         _error(3360, "failure to converge for element " + strofreal(i))
        if (rc)
         _error(3498, "invalid search range for element " + strofreal(i))
        res[i] = resi
    }
    return(res)

}

real scalar _mm_finvert_brent_f(
    real scalar x,
    real scalar y,
    pointer(real scalar function) scalar f,
    | opt)
{
    if (args()<4)
        return( (*f)(x) - y )
    else
        return( (*f)(x, opt) - y )
}

end


*! {smcl}
*! {marker mm_integrate}{bf:mm_integrate.mata}{asis}
*! version 1.0.0  Ben Jann  26jan2014

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

// Simpson's rule (quadratic interpolation)
real colvector mm_integrate_sr(
    pointer(real scalar function) scalar f, // function to be integrated
    real scalar a,      // lower limit
    real scalar b,      // upper limit
    real scalar n,      // integration intervals
    | real scalar el,   // element wise evaluation
    `opts')             // additional args to pass on to f
{
    transmorphic    fs            // setup for f
    real scalar     d, i
    real colvector  x, w
    
    // check options
    if (missing(n))       _error(3300, "n: missing not allowed")
    if ((n<2) | mod(n,2)) _error(3300, "n: must be a positive multiple of two")
    if (missing(a))       _error(3300, "a: missing not allowed")
    if (missing(b))       _error(3300, "b: missing not allowed")
    if (a>=b)             _error(3300, "a must be strictly smaller than b")
    
    // prepare function call
    fs = mm_callf_setup(f, args()-5, `opts')
    
    // integration setup
    d = (b-a)/n
    x = rangen(a, b, n+1)
    w = 1 \ colshape((J(n/2,1,4), J(n/2,1,2)), 1)
    w[n+1] = 1
    
    // pass column vector
    if (el==0 | args()<5) {
        return(d/3 * quadcolsum(mm_callf(fs, x) :* w))
    }
    // element wise evaluation
    for (i=1; i<=(n+1); i++) {
        x[i] = mm_callf(fs, x[i])
    }
    return(d/3 * quadcolsum(x :* w))
}

// Simpson's 3/8 rule (cubic interpolation)
real colvector mm_integrate_38(
    pointer(real scalar function) scalar f, // function to be integrated
    real scalar a,      // lower limit
    real scalar b,      // upper limit
    real scalar n,      // integration intervals
    | real scalar el,   // element wise evaluation
    `opts')             // additional args to pass on to f
{
    transmorphic    fs            // setup for f
    real scalar     d, i
    real colvector  x, w
    
    // check options
    if (missing(n))       _error(3300, "n: missing not allowed")
    if ((n<3) | mod(n,3)) _error(3300, "n: must be a positive multiple of three")
    if (missing(a))       _error(3300, "a: missing not allowed")
    if (missing(b))       _error(3300, "b: missing not allowed")
    if (a>=b)             _error(3300, "a must be strictly smaller than b")
    
    // prepare function call
    fs = mm_callf_setup(f, args()-5, `opts')
    
    // integration setup
    d = (b-a)/n
    x = rangen(a, b, n+1)
    w = 1 \ colshape((J(n/3,1,3), J(n/3,1,3), J(n/3,1,2)), 1)
    w[n+1] = 1
    
    // pass column vector
    if (el==0 | args()<5) {
        return(3*d/8 * quadcolsum(mm_callf(fs, x) :* w))
    }
    // element wise evaluation
    for (i=1; i<=(n+1); i++) {
        x[i] = mm_callf(fs, x[i])
    }
    return(3*d/8 * quadcolsum(x :* w))
}

end


*! {smcl}
*! {marker mm_ipolate}{bf:mm_ipolate.mata}{asis}
*! version 2.0.0  08jul2020  Ben Jann
version 9.2
mata:

real colvector mm_ipolate(real colvector x0, real colvector y0,
    real colvector x1, | real scalar outer)
{
    real colvector p0, p1, y1
    
    if (args()<4) outer = 0
    if (rows(y0)!=rows(x0)) _error(3200)
    p0 = order(x0,1)
    p1 = order(x1,1)
    y1 = J(rows(x1),1,.)
    y1[p1] = _mm_ipolate(x0[p0], y0[p0], x1[p1], outer)
    return(y1)
}

real colvector _mm_ipolate(real colvector x0, real colvector y0,
    real colvector x1, | real scalar outer)
{   // x0 and x1 assumed sorted
    real scalar    n, i, ai, bi
    real colvector a, b, x, y

    if (args()<4) outer = 0
    n = rows(x0)
    if (rows(y0)!=n) _error(3200)
    if (n<=1) return(mm_fastipolate(x0, y0, x1, outer))
    // collapse data
    a = x0:!=(x0[n]\x0[|1\n-1|]) // tag first obs within panel
    a[1] = 1                     // (should x0 be constant)
    a = select(1::n, a)          // get index of first obs within panel
    if ((i=rows(a))==n) return(mm_fastipolate(x0, y0, x1, outer)) // no ties
    b = x0:!=(x0[|2\.|]\x0[1])   // tag last obs within panel
    b[n] = 1                     // (should x0 be constant)
    b = select(1::n, b)          // get index of last obs within panel
    x = y = J(i,1,.)
    for (;i;i--) {
        ai = a[i]; bi = b[i]
        x[i] = x0[ai]
        if (ai==bi) {
            y[i] = y0[ai]
            continue
        }
        y[i] = mean(y0[|ai\bi|],1)
    }
    return(mm_fastipolate(x, y, x1, outer))
}

real colvector mm_fastipolate(real colvector x0, real colvector y0,
    real colvector x1, | real scalar outer)
{   // x0 and x1 assumed sorted; x0 assumed unique
    real scalar    n, i, j
    real colvector y1

    if (args()<4) outer = 0
    j = rows(x0); n = i = rows(x1)
    if (rows(y0)!=j) _error(3200)
    if (j<1) return(J(n, 1, .))
    y1 = J(n,1,.)
    // handle values of x1 > max(x0)
    for (;i;i--) {
        if (x1[i]<=x0[j]) break
    }
    if (outer) {
        if (i<n) y1[|i+1\.|] = J(n-i, 1, y0[j])
    }
    // handle values of x1 within minmax(x0)
    for (;i;i--) {
        for (;j;j--) {
            if (x0[j]==x1[i]) {
                y1[i] = y0[j]
                break
            }
            if (x0[j]<x1[i]) {
                y1[i] = y0[j] + (y0[j+1]-y0[j]) * (x1[i]-x0[j])/(x0[j+1]-x0[j])
                break
            }
        }
        if (!j) break
    }
    // handle values of x1 < min(x0)
    if (outer) {
        if (i) y1[|1\i|] = J(i, 1, y0[1])
    }
    return(y1)
}

end


*! {smcl}
*! {marker mm_polint}{bf:mm_polint.mata}{asis}
*! version 1.0.0, Ben Jann, 12jul2006
version 9.2
mata:

real vector mm_polint(
 real vector xa,
 real vector ya,
 real vector x,
 | real scalar degree) // degree
{
	real scalar  n, m, i, j, a, b
	real vector  y

	if ( args()<4 ) degree = 1
	n = length(xa)
	if (length(ya)!=n) _error(3200)
	if ( degree<1 ) _error(3498, "degree must be 1 or larger")
	if ( degree>=n ) _error(3498, "degree must be smaller than length of data")
	if ( degree!=trunc(degree) ) _error(3498, "degree must be integer")

	m = degree + 1
	y = J(rows(x), cols(x), .)
	j = 0
	for (i=1; i<=length(x); i++) {
		mm_hunt(xa, x[i], j)
		a = min((max((trunc(j-(m-1)/2+.5), 1)), n+1-m))
		b = a+m-1
		y[i] = _mm_polint(xa[|a \ b|], ya[|a \ b|], x[i])
	}
	return(y)
}

real scalar _mm_polint(
// translation of -polint- from Press et al. (1992), Numerical
// Recipes in C, p. 109-110, http://www.numerical-recipes.com
 real vector xa,
 real vector ya,
 real scalar x)
{
	real scalar  i, m, n, ns
	real scalar  y, den, dif, dift, ho, hp, w
	real vector  c, d

	n = length(xa)
	if ( length(ya) != n ) _error(3200)

	ns = 1
	dif = abs(x-xa[1])
	c = J(1, n, .)
	d = J(1, n, .)
	for (i=1; i<=n; i++) {
		if ( (dift=abs(x-xa[i])) < dif) {
			ns  = i
			dif = dift
		}
		c[i] = ya[i]
		d[i] = ya[i]
	}
	y = ya[ns--]
	for (m=1; m<n; m++) {
		for (i=1; i<=n-m; i++) {
			ho = xa[i] - x
			hp = xa[i+m] - x
			w = c[i+1] - d[i]
			if ( (den=ho-hp) == 0.0)
			 _error(3498, "invalid input data") // This error can occur
			 // only if two input xa's are (to within roundoff) identical.
			den = w / den
			d[i] = hp * den
			c[i] = ho * den
		}
		if ( 2*ns < (n-m) ) y = y + c[ns+1]
		else y = y + d[ns--]
	}
	return(y)
}

end


*! {smcl}
*! {marker mm_sqrt}{bf:mm_sqrt.mata}{asis}
*! version 1.0.0  29may2017  Christopher Baum and Ben Jann
version 9.2
mata:

real matrix mm_sqrt(real matrix A)
{
    real colvector  L
    real matrix     X
    pragma unset    L
    pragma unset    X
    
    if (isfleeting(A)) _symeigensystem(A, X, L)
    else                symeigensystem(A, X, L)
    return(X * diag(sqrt(L)) * X')
}

end


*! {smcl}
*! {marker mm_plot}{bf:mm_plot.mata}{asis}
*! version 1.0.0, Ben Jann, 03aug2006
version 9.2
mata:

void mm_plot(
 real matrix X,
 | string scalar type,
   string scalar opts)
{
	real scalar rc

	rc = _mm_plot(X, type, opts)
	if (rc) _error(rc)
}

real scalar _mm_plot(
 real matrix X,
 | string scalar type,
   string scalar opts)
{
	real scalar   n, N, k, j, rc, updata
	string scalar vars, plottype
	real vector   vid

	updata = st_updata()
	plottype = type
	if (plottype=="") plottype = "scatter"
	n = rows(X)
	k = cols(X)
	N = st_nobs()
	if (N<n & k>0) st_addobs(n-N)
	vid = J(1, k, .)
	for (j=1; j<=k; j++) {
		st_store((1,n), vid[j]=st_addvar("double", st_tempname()), X[,j])
		st_varlabel(vid[j], "Column " + strofreal(j))
		vars = vars + " " + st_varname(vid[j])
	}
	st_updata(updata)
	rc = _stata("twoway " + plottype + vars + ", " + opts)
	updata = st_updata() // twoway might change data-have-changed flag
	if (N<n & k>0) st_dropobsin((N+1,n))
	st_dropvar(vid)
	st_updata(updata)
	return(rc)
}

end


*! {smcl}
*! {marker mm_group}{bf:mm_group.mata}{asis}
*! version 1.0.0  09jul2020  Ben Jann
version 9.2
mata:

real colvector mm_group(transmorphic matrix X, | real rowvector idx)
{
    real scalar    r, c
    real colvector p
    
    r = rows(X); c = cols(X)
    if (args()<2) {
        if (r<=1 | c==0) return(J(r,1,1))
        p = order(X, 1..c)
        p[p] = _mm_group(X[p,])
    }
    else {
        if (length(idx)>c) _error(3200)
        if (r<=1 | c==0) return(J(r,1,1))
        p = order(X, idx)
        p[p] = _mm_group(X[p, abs(idx)])
    }
    return(p)
}

real colvector _mm_group(transmorphic matrix X)
{
    real scalar r, c
    real colvector p
    
    r = rows(X); c = cols(X)
    if (r<=1 | c==0) return(J(r,1,1))
    if (c==1) p = (X :!= (X[r] \ X[|1\r-1|]))
    else      p = (rowsum(X :!= (X[r,] \ X[|1,1 \ r-1,.|])) :!= 0)
    return(mm_colrunsum(p))
}

end




*! {smcl}
*! {marker mm_panels}{bf:mm_panels.mata}{asis}
*! version 1.0.3, Ben Jann, 03may2007
version 9.2
mata:

void function mm_panels(
 transmorphic vector X1,       // upper level ID variable (e.g. strata)
 transmorphic matrix info1,    // will be replaced (unless X1 and X2 missing)
 | transmorphic vector X2,     // lower level ID variable (e.g. clusters)
   transmorphic matrix info2)  // will be replaced (unless X2 missing)
{
	real scalar i, j, nc, b, e, b2, e2

	if (length(X1)>0 & X1!=.) info1 = _mm_panels(X1)
	if (args()<3) return
	if (length(X2)==0 | X2==.) {
		if (length(X1)>0 & X1!=.) info1 = info1, info1
		return
	}
	if (length(X1)==0 | X1==.) {
		info2 = _mm_panels(X2)
		info1 = colsum(info2), rows(info2)
		return
	}
	if (length(X1)!=length(X2)) _error(3200)
	info1 = info1, J(rows(info1),1,.)
	e = 0
	for (i=1; i<=rows(info1); i++) {
		nc = 1
		b = e + 2
		e = e + info1[i,1]
		for (j=b; j<=e; j++) {
			if (X2[j]!=X2[j-1]) nc++
		}
		info1[i,2] = nc
	}
	if (args()<4) return
	info2 = J(colsum(info1[.,2]), 1, .)
	e = e2 = 0
	for (i=1; i<=rows(info1); i++) {
		b = e + 1
		e = e + info1[i,1]
		b2 = e2 + 1
		e2 = e2 + info1[i,2]
		info2[|b2 \ e2|] = _mm_panels(X2[|b \ e|], info1[i,2])
	}
}

real colvector _mm_panels(transmorphic vector X, | real scalar np)
{
	real scalar i, j, r
	real colvector res

	r = length(X)
	if (r<1) return(J(0,1,.))
	if (args()<2) np = r
	res = J(np, 1, 1)
	j = 1
	for (i=2; i<=r; i++) {
		if (X[i]==X[i-1]) res[j] = res[j] + 1
		else j++
	}
	if (j==r) return(res)
	return(res[|1 \ j|])
}

// /*Old (slow) version of _mm_panels()*/
//real colvector _mm_panels(transmorphic vector X, | real scalar np)
//{
//	real scalar i, j, n
//	real colvector res
//
//	if (args()<2) np = mm_npanels(X)
//	if (length(X)==0) return(J(0,1,.))
//	res = J(np, 1, .)
//	n = j = 1
//	for (i=2; i<=length(X); i++) {
//		if (X[i]!=X[i-1]) {
//			res[j++] = n
//			n = 1
//		}
//		else n++
//	}
//	res[j] = n
//	return(res)
//}

real scalar mm_npanels(vector X)
{
	real scalar i, np

	if (length(X)==0) return(0)
	np = 1
	for (i=2; i<=length(X); i++) {
		if (X[i]!=X[i-1]) np++
	}
	return(np)
}

end


*! {smcl}
*! {marker mm_nunique}{bf:mm_nunique.mata}{asis}
*! version 1.1.1  09oct2020  Ben Jann
version 9.2
mata:

real scalar mm_nunique(transmorphic vector X)
{
    return(sum(mm_unique_tag(X)))
}

transmorphic vector mm_unique(transmorphic vector X, | real scalar order)
{
    if (rows(X)*cols(X)==0)  return(X)
    if (order==1 | order==2) return(select(X, mm_unique_tag(X, order)))
    // sorted
    if (cols(X)==1)    return(sort(select(X, mm_unique_tag(X)),1))
    if (!iscomplex(X)) return(sort(select(X, mm_unique_tag(X))',1)')
    return(transposeonly(sort(transposeonly(select(X, mm_unique_tag(X))),1)))
}

/*
real vector mm_unique_idx(transmorphic vector X, | real scalar order)
{
    if (rows(X)*cols(X)==0) return(J(rows(X),cols(X),.))
    if (cols(X)==1) return(select(1::rows(X), mm_unique_tag(X, order)))
                    return(select(1..cols(X), mm_unique_tag(X, order)))
}
*/

real vector mm_unique_tag(transmorphic vector X, | real scalar order)
{
    real vector p
    
    if (rows(X)*cols(X)==0) return(J(rows(X),cols(X),.))
    if (order==1) {
        // use stable sort order and tag first occurrence
        if (cols(X)==1)        p = mm_order(X,1,1)
        else if (iscomplex(X)) p = mm_order(transposeonly(X),1,1)'
        else                   p = mm_order(X',1,1)'
        p[p] = _mm_unique_tag(X[p])
    }
    else if (order==2) {
        // use stable sort order and tag last occurrence
        if (cols(X)==1)        p = mm_order(X,1,1)
        else if (iscomplex(X)) p = mm_order(transposeonly(X),1,1)'
        else                   p = mm_order(X',1,1)'
        p[p] = _mm_unique_tag(X[p], 1)
    }
    else {
        // use non-stable sort order (tag random element within ties)
        if (cols(X)==1)        p = order(X,1)
        else if (iscomplex(X)) p = order(transposeonly(X),1)'
        else                   p = order(X',1)'
        p[p] = _mm_unique_tag(X[p])
    }
    return(p)
}

real scalar _mm_nunique(transmorphic vector X)
{
    return(sum(_mm_unique_tag(X)))
}

transmorphic vector _mm_unique(transmorphic vector X)
{
    if (rows(X)*cols(X)==0) return(X)
    return(select(X, _mm_unique_tag(X)))
}

/*
real vector _mm_unique_idx(transmorphic vector X, | real scalar last)
{
    if (args()<2) last = 0
    if (rows(X)*cols(X)==0) return(J(rows(X),cols(X),.))
    if (cols(X)==1) return(select(1::rows(X), _mm_unique_tag(X, last)))
                    return(select(1..cols(X), _mm_unique_tag(X, last)))
}
*/

real vector _mm_unique_tag(transmorphic vector X, | real scalar last)
{
    real scalar r, c
    real vector p

    if (args()<2) last = 0
    r = rows(X); c = cols(X)
    if (r*c==0)      return(J(r,c,.))
    if (r==1 & c==1) return(J(1,1,1))
    if (last) {
        if (c==1) {
            p = (X :!= (X[|2\.|] \ X[1]))
            p[r] = 1 // should X be constant
            return(p)
        }
        p = (X :!= (X[|2\.|] , X[1]))
        p[c] = 1 // should X be constant
        return(p)
    }
    if (c==1) {
        p = (X :!= (X[r] \ X[|1\r-1|]))
        p[1] = 1 // should X be constant
        return(p)
    }
    p = (X :!= (X[c] , X[|1\c-1|]))
    p[1] = 1 // should X be constant
    return(p)
}

real scalar mm_nuniqrows(transmorphic matrix X)
{
    return(sum(mm_uniqrows_tag(X)))
}

transmorphic matrix mm_uniqrows(transmorphic matrix X, | real scalar order)
{
    if (rows(X)==0) return(X)
    if (order==1 | order==2) return(select(X, mm_uniqrows_tag(X, order)))
    return(mm_sort(select(X, mm_uniqrows_tag(X))))
}

/*
real vector mm_uniqrows_idx(transmorphic matrix X, | real scalar order)
{
    if (rows(X)==0) return(J(0,1,.))
    return(select(1::rows(X), mm_uniqrows_tag(X, order)))
}
*/

real vector mm_uniqrows_tag(transmorphic matrix X, | real scalar order)
{
    real vector p
    
    if (rows(X)==0) return(J(0,1,.))
    if (order==1) {
        p = mm_order(X,.,1)
        p[p] = _mm_uniqrows_tag(X[p,])
    }
    else if (order==2) {
        p = mm_order(X,.,1)
        p[p] = _mm_uniqrows_tag(X[p,], 1)
    }
    else {
        p = mm_order(X,.,0)
        p[p] = _mm_uniqrows_tag(X[p,])
    }
    return(p)
}

real scalar _mm_nuniqrows(transmorphic matrix X)
{
    return(sum(_mm_uniqrows_tag(X)))
}

transmorphic matrix _mm_uniqrows(transmorphic matrix X)
{
    if (rows(X)==0) return(X)
    return(select(X, _mm_uniqrows_tag(X)))
}

/*
real vector _mm_uniqrows_idx(transmorphic matrix X, | real scalar last)
{
    if (args()<2) last = 0
    if (rows(X)==0) return(J(0,1,.))
    return(select(1::rows(X), _mm_uniqrows_tag(X, last)))
}
*/

real vector _mm_uniqrows_tag(transmorphic matrix X, | real scalar last)
{
    real scalar r, c
    real vector p

    if (args()<2) last = 0
    r = rows(X); c = cols(X)
    if (r==0) return(J(0,1,.))
    if (r==1) return(J(1,1,1))
    if (last) {
        if (c==0) p = J(r,1,0)
        else      p = (rowsum(X :!= (X[|2,1 \ .,.|] \ X[1,])):!=0)
        p[r] = 1 // should X be constant
        return(p)
    }
    if (c==0) p = J(r,1,0)
    else      p = (rowsum(X :!= (X[r,] \ X[|1,1 \ r-1,.|])):!=0)
    p[1] = 1 // should X be constant
    return(p)
}

end



*! {smcl}
*! {marker mm_diff}{bf:mm_diff.mata}{asis}
*! version 1.0.0  Ben Jann  07jul2020
version 9.2
mata:

numeric vector mm_diff(numeric vector x, | real scalar lag0)
{
    real scalar n, lag
    
    if (lag0>=.) lag =  1
    else         lag = abs(trunc(lag0))
    n = length(x)
    if (n<=lag) {
        if (cols(x)!=1) return(J(1,0,missingof(x)))
        return(J(0,1,missingof(x)))
    }
    return(x[|1+lag \ .|] :- x[|1 \ n-lag|])
}

numeric matrix mm_rowdiff(numeric matrix x, | real scalar lag0)
{
    real scalar n, lag
    
    if (lag0>=.) lag =  1
    else         lag = abs(trunc(lag0))
    n = cols(x)
    if (n<=lag)     return(J(rows(x),0,missingof(x)))
    if (rows(x)==0) return(J(0, n-lag, missingof(x)))
    return(x[|1,1+lag \ .,.|] :- x[|1,1 \ .,n-lag|])
}

numeric matrix mm_coldiff(numeric matrix x, | real scalar lag0)
{
    real scalar n, lag
    
    if (lag0>=.) lag =  1
    else         lag = abs(trunc(lag0))
    n = rows(x)
    if (n<=lag)     return(J(0,cols(x),missingof(x)))
    if (cols(x)==0) return(J(n-lag, 0, missingof(x)))
    return(x[|1+lag,1 \ .,.|] :- x[|1,1 \ n-lag,.|])
}

end


*! {smcl}
*! {marker mm_isconstant}{bf:mm_isconstant.mata}{asis}
*! version 1.0.2, Ben Jann, 16jul2020
version 9.0
mata:

real scalar mm_isconstant(transmorphic matrix X)
{
    if (length(X)<2) return(1)
    return(allof(X, X[1,1]))
}

real scalar mm_issorted(transmorphic vector x, | real scalar descending)
{
    if (ispointer(x)) _error(3250)
    if (length(x)<=1) return(1)
    if (args()<2) descending = 0
    if (descending) {
        if (isstring(x)) {
            if (cols(x)==1) return(all(x[|1 \ length(x)-1|] :>= x[|2 \ .|]))
                            return(all(x[|1 \ length(x)-1|] :>= x[|2 \ .|]))
        }
        if (cols(x)==1) return(all((.z \ x[|1 \ length(x)-1|]) :>= x))
                        return(all((.z, x[|1 \ length(x)-1|]) :>= x))
    }
    if (isstring(x)) {
        if (cols(x)==1) return(all(x[|1 \ length(x)-1|] :<= x[|2 \ .|]))
                        return(all(x[|1 \ length(x)-1|] :<= x[|2 \ .|]))
    }
    if (cols(x)==1) return(all(x :<= (x[|2 \ .|] \ .z)))
                    return(all(x :<= (x[|2 \ .|], .z)))
}

end


*! {smcl}
*! {marker mm_colrunsum}{bf:mm_colrunsum.mata}{asis}
*! version 1.0.9  09jul2020  Ben Jann
version 9.0
mata:

numeric matrix mm_colrunsum(numeric matrix A, | real scalar mis, real scalar qd)
{
    numeric matrix B

    if (args()<2) mis = 0
    if (stataversion()>=1000) {
        if (args()<3) qd = 0
        return(_mm_colrunsum_goto10(A, mis, qd))
    }
    if (isfleeting(A)) {
        _mm_colrunsum(A, mis)
        return(A)
    }
    _mm_colrunsum(B=A, mis)
    return(B)
}

void _mm_colrunsum(numeric matrix Z, real scalar mis)
{
    real scalar i

    if (mis==0) _editmissing(Z, 0)
    for (i=2; i<=rows(Z); i++) Z[i,] = Z[i-1,] + Z[i,]

//    // running sum using mean-update formula; see Gould 2006, SJ 6(4)
//    if (mis==0) _editmissing(Z, 0)
//    if (rows(Z)<2) return
//    for (i=2; i<=rows(Z); i++) Z[i,] = Z[i-1,] + (Z[i,]-Z[i-1,])/i
//    Z = Z :* (1::rows(Z))

}

numeric matrix _mm_colrunsum_goto10(numeric matrix Z, real scalar mis, real scalar qd)
{
    if (qd) return(_mm_quadcolrunsum10(Z, mis))
    return(_mm_colrunsum10(Z, mis))
}


end


*! {smcl}
*! {marker mm_linbin}{bf:mm_linbin.mata}{asis}
*! version 1.0.9  12aug2020  Ben Jann
version 9.0
mata:

real colvector mm_linbin(
    real colvector x,  // data points
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real colvector  p

    p = mm_order(x, ., 1) // stable sort order
    if (rows(w)==1) return(__mm_linbin(x[p], g)*w)
                    return(__mm_linbin_w(x[p], w[p], g))
}

real colvector _mm_linbin(
    real colvector x,  // data points (assumed sorted)
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    if (rows(w)==1) return(__mm_linbin(x, g)*w)
                    return(__mm_linbin_w(x, w, g))
}

real colvector __mm_linbin(
    real colvector x,  // data points (assumed sorted)
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     i, j, g1, g0, i1, i0
    real colvector  c, xx

    j = rows(g)
    i = rows(x)
    if (i<1) return(J(j, 1, 0))
    if (j<2) return(i)
    c = J(j, 1, 0)
    // data above grid
    g0 = g[j]
    i1 = i
    for (;i;i--) {
        if (x[i]<g0) break
    }
    i0 = i+1
    c[j] = (i1-i0+1)
    // data within grid range
    g1 = g0
    for (j--; j; j--) {
        g1 = g0
        g0 = g[j]
        i1 = i
        for (;i;i--) {
            if (x[i]<g0) break
        }
        i0 = i+1
        if (i1<i0) continue // no data in current bin
        xx = x[|i0\i1|]
        c[j] =            sum(g1 :- xx) / (g1 - g0)
        c[j+1] = c[j+1] + sum(xx :- g0) / (g1 - g0)
    }
    // data below grid
    if (i) c[1] = c[1] + i
    return(c)
}

real colvector __mm_linbin_w(
    real colvector x,  // data points (assumed sorted)
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     i, j, g1, g0, i1, i0
    real colvector  c, xx, ww

    j = rows(g)
    i = rows(x)
    if (rows(w)!=i) _error(3200)
    if (i<1) return(J(j, 1, 0))
    if (j<2) return(sum(w))
    c = J(j, 1, 0)
    // data above grid
    g0 = g[j]
    i1 = i
    for (;i;i--) {
        if (x[i]<g0) break
    }
    i0 = i+1
    if (i0<=i1) c[j] = sum(w[|i0 \ i1|])
    // data within grid range
    g1 = g0
    for (j--; j; j--) {
        g1 = g0
        g0 = g[j]
        i1 = i
        for (;i;i--) {
            if (x[i]<g0) break
        }
        i0 = i+1
        if (i1<i0) continue // no data in current bin
        xx = x[|i0\i1|]
        ww = w[|i0\i1|]
        c[j] =            sum(ww:*(g1 :- xx)) / (g1 - g0)
        c[j+1] = c[j+1] + sum(ww:*(xx :- g0)) / (g1 - g0)
    }
    // data below grid
    if (i) c[1] = c[1] + sum(w[|1 \ i|])
    return(c)
}

real colvector mm_fastlinbin(
    real colvector x,  // data points
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     g1, i, j, d, z, ng
    real colvector  c

    ng = rows(g)
    i  = rows(x)
    if (ng<1) _error(3200)
    if (i<1)  return(J(ng, 1, 0))
    if (ng<2) return(mm_nobs(x,w))
    if (rows(w)!=1) return(_mm_fastlinbin_w(x, w, g))
    g1 = g[1]
    d = (g[ng] - g1)/(ng-1)
    c = J(ng, 1, 0)
    for (; i; i--) {
        z = (x[i] - g1) / d
        if (z<=0) {
            c[1] = c[1] +  1
            continue
        }
        if (z>=(ng-1)) {
            c[ng] = c[ng] +  1
            continue
        }
        j = trunc(z) + 1
        c[j] = c[j] + (j-z)
        j++
        c[j] = c[j] + (2-j+z)
    }
    return(c*w)
}

real colvector _mm_fastlinbin_w(
    real colvector x,  // data points
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     g1, i, j, d, z, ng
    real colvector  c

    ng = rows(g)
    i  = rows(x)
    if (rows(w)!=i) _error(3200)
    g1 = g[1]
    d = (g[ng] - g1)/(ng-1)
    c = J(ng, 1, 0)
    for (; i; i--) {
        z = (x[i] - g1) / d
        if (z<=0) {
            c[1] = c[1] +  w[i]
            continue
        }
        if (z>=(ng-1)) {
            c[ng] = c[ng] +  w[i]
            continue
        }
        j = trunc(z) + 1
        c[j] = c[j] + w[i] * (j-z)
        j++
        c[j] = c[j] + w[i] * (2-j+z)
    }
    return(c)
}

end


*! {smcl}
*! {marker mm_exactbin}{bf:mm_exactbin.mata}{asis}
*! version 1.0.6  12aug2020  Ben Jann
version 9.0
mata:

real colvector mm_exactbin(
    real colvector x,   // data
    real colvector w,   // weights
    real colvector g,   // grid (assumed sorted)
    | real scalar dir,  // 0 left inclusive (default),  else right inclusive
      real scalar include) // include data outside grid in first/last bin
{
    real colvector  p

    if (args()<4) dir = 0
    if (args()<5) include = 0
    p = mm_order(x, ., 1) // stable sort order
    if (rows(w)==1) return(__mm_exactbin(x[p], g, dir, include)*w)
                    return(__mm_exactbin_w(x[p], w[p], g, dir, include))
}

real colvector _mm_exactbin(
    real colvector x,   // data (assumed sorted)
    real colvector w,   // weights
    real colvector g,   // grid (assumed sorted)
    | real scalar dir,  // 0 left inclusive (default),  else right inclusive
      real scalar include) // include data outside grid in first/last bin
{
    if (args()<4) dir = 0
    if (args()<5) include = 0
    if (rows(w)==1) return(__mm_exactbin(x, g, dir, include)*w)
                    return(__mm_exactbin_w(x, w, g, dir, include))
}

real colvector __mm_exactbin(
    real colvector x,    // data (assumed sorted)
    real colvector g,    // grid (assumed sorted)
    real scalar dir,     // 0 left inclusive (default),  else right inclusive
    real scalar include) // include data outside grid in first/last bin
{
    real scalar     i, i1, j, gj
    real colvector  c

    j = rows(g)
    if (j<2) _error(3200)
    i = rows(x)
    if (i<1) return(J(j-1, 1, 0))
    if (include==0) {
        if (x[1]<g[1] | x[i]>g[j]) _error("data out of grid range")
    }
    if (j<3) return(i)
    j--
    c = J(j, 1, 0)
    if (dir==0) {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = i1 - i
        }
    }
    else {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<=gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = i1 - i
        }
    }
    if (i) c[1] = c[1] + i // data below grid
    return(c)
}

real colvector __mm_exactbin_w(
    real colvector x,    // data (assumed sorted)
    real colvector w,    // weights
    real colvector g,    // grid (assumed sorted)
    real scalar dir,     // 0 left inclusive (default),  else right inclusive
    real scalar include) // include data outside grid in first/last bin
{
    real scalar     i, i1, j, gj
    real colvector  c
    
    j = rows(g)
    if (j<2) _error(3200)
    i = rows(x)
    if (i<1) return(J(j-1, 1, 0))
    if (include==0) {
        if (x[1]<g[1] | x[i]>g[j]) _error("data out of grid range")
    }
    if (j<3) return(sum(w))
    j--
    c = J(j, 1, 0)
    if (dir==0) {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = sum(w[|i+1 \ i1|])
        }
    }
    else {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<=gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = sum(w[|i+1 \ i1|])
        }
    }
    if (i) c[1] = c[1] + sum(w[|1 \ i|]) // data below grid
    return(c)
}

real colvector mm_fastexactbin(
    real colvector x,      // data points
    real colvector w,      // weights
    real colvector g,      // grid (assumed sorted)
    | real scalar dir,     // 0 left inclusive (default),  else right inclusive
      real scalar include) // include data outside grid in first/last bin
{
    real scalar     g1, i, j, n, d, xi
    real colvector  c
    real rowvector  minmax
    pointer scalar  wi

    if (args()<4) dir = 0
    if (args()<5) include = 0
    n = rows(g) - 1
    i = rows(x)
    if (i<1) return(J(n, 1, 0))
    if (include==0) {
        minmax = minmax(x)
        if (minmax[1]<g[1] | minmax[2]>g[n+1]) _error("data out of grid range")
    }
    if (n<2) return(mm_nobs(x,w))
    wi = (rows(w)==1 ? &1 : (rows(w)!=i ? _error(3200) : &i))
    g1 = g[1]
    d = (g[n+1] - g1)/n
    c = J(n, 1, 0)
    if (dir==0) {
        for (; i; i--) {
            xi = x[i]
            j = floor((xi - g1) / d) + 1
            // data below grid
            if (j<1) {
                c[1] = c[1] +  w[*wi]
                continue
            }
            // data above grid
            if (j>n) {
                c[n] = c[n] +  w[*wi]
                continue
            }
            // roundoff error 1: j erroneously rounded up
            if (xi<g[j]) {
                // not sure whether this can happen at all...
                if (j>1) c[j-1] = c[j-1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // roundoff error 2: j erroneously rounded down
            if (xi>=g[j+1]) {
                if (j<n) c[j+1] = c[j+1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // no roundoff error
            c[j] = c[j] + w[*wi]
        }
    }
    else {
        for (; i; i--) {
            xi = x[i]
            j = ceil((xi - g1) / d)
            // data below grid
            if (j<1) {
                c[1] = c[1] +  w[*wi]
                continue
            }
            // data above grid
            if (j>n) {
                c[n] = c[n] +  w[*wi]
                continue
            }
            // roundoff error 1: j erroneously rounded up
            if (xi<=g[j]) {
                if (j>1) c[j-1] = c[j-1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // roundoff error 2: j erroneously rounded down
            if (xi>g[j+1]) {
                // not sure whether this can happen at all...
                if (j<n) c[j+1] = c[j+1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // no roundoff error
            c[j] = c[j] + w[*wi]
        }
    }
    return(c)
}

end


*! {smcl}
*! {marker mm_makegrid}{bf:mm_makegrid.mata}{asis}
*! version 1.0.4, Ben Jann, 26jun2006
version 9.0
mata:

real colvector mm_makegrid(
 real colvector x, // data points
 | real scalar M,  // number of bins (default: 401)
   real scalar e,  // extend range by e on each side
   real scalar min,  // min
   real scalar max)  // max
{
	real scalar a, b

	a = (min<. ? min : min(x) - (e<. ? e : 0))
	b = (max<. ? max : max(x) + (e<. ? e : 0))
	return( rangen(a, b, (M<. ? M : 512)) )
}

end


*! {smcl}
*! {marker mm_seq}{bf:mm_seq.mata}{asis}
*! version 1.0.0  Ben Jann  04aug2020

version 9.2
mata:

real colvector mm_seq(real scalar a, real scalar b, real scalar d)
{
    if (missing((a,b,d))) return(J(0,1,.))
    if (a<=b)             return( a :+ (0::(b-a)/abs(d)) * abs(d) )
                          return( a :- (0::(a-b)/abs(d)) * abs(d) )
}

end


*! {smcl}
*! {marker mm_cut}{bf:mm_cut.mata}{asis}
*! version 1.0.0, Ben Jann, 16mar2006
version 9.1
mata:

real colvector mm_cut(real colvector x, real vector at, | real scalar sorted)
{
	real scalar i, j
	real colvector result, p

	result = J(rows(x),1,.)
	j = length(at)
	if (args()<3) sorted = 0
	if (sorted) {
		for (i=rows(x); i>0; i--) {
			if (x[i]>=.) continue
			for (; j>0; j--) {
				if (at[j]<=x[i]) break
			}
			if (j>0) result[i,1] = at[j]
		}
		return(result)
	}
	p = order(x,1)
	for (i=rows(x); i>0; i--) {
		if (x[p[i]]>=.) continue
		for (; j>0; j--) {
			if (at[j]<=x[p[i]]) break
		}
		if (j>0) result[p[i],1] = at[j]
	}
	return(result)
}
end


*! {smcl}
*! {marker mm_posof}{bf:mm_posof.mata}{asis}
*! version 1.0.0, Ben Jann, 24mar2006
version 9.0
mata:

real scalar mm_posof(transmorphic vector haystack, transmorphic scalar needle)
{
	real scalar i

	for (i=1; i<=length(haystack); i++) {
		if (haystack[i]==needle) return(i)
	}
	return(0)
}
end


*! {smcl}
*! {marker mm_which}{bf:mm_which.mata}{asis}
*! version 1.0.2, Ben Jann, 17apr2007
version 9.2
mata:

real matrix mm_which(real vector I)
{
	if (cols(I)!=1) return(select(1..cols(I), I))
	else return(select(1::rows(I), I))
}

end


*! {smcl}
*! {marker mm_locate}{bf:mm_locate.mata}{asis}
*! version 1.0.0, Ben Jann, 07jul2006

version 9.2

mata:

void mm_locate(
// translation of -locate- from Press et al. (1992), Numerical
// Recipes in C, p. 117, http://www.numerical-recipes.com
// Description: "Given an array xx[1..n], and given a value x,
// returns a value j such that x is between xx[j] and xx[j+1].
// xx must be monotonic, either increasing or decreasing. j=0
// or j=n is returned to indicate that x is out of range."
 numeric vector xx,
 numeric scalar x,
 real scalar j)
{
	real scalar  ju, jm, jl
	real scalar  ascnd
	real scalar  n

	n = length(xx)
	jl = 0
	ju = n + 1
	ascnd = (xx[n] >= xx[1])
	while (ju-jl > 1) {
		jm = trunc((ju+jl)/2)
		if (x >= xx[jm] == ascnd)
		  jl = jm
		else
		  ju = jm
	}
	if (x == xx[1])     j = 1
	else if(x == xx[n]) j = n - 1
	else                j = jl
}

void mm_hunt(
// translation of -hunt- from Press et al. (1992), Numerical
// Recipes in C, p. 118-9, http://www.numerical-recipes.com
// Description: "Given an array xx[1..n], and given a value x,
// returns a value jlo such that x is between xx[jlo] and
// xx[jlo+1]. xx[1..n] must be monotonic, either increasing or
// decreasing. jlo=0 or jlo=n is returned to indicate that x is
// out of range. jlo on input is taken as the initial guess for
// jlo on output.
 numeric vector xx,
 numeric scalar x,
 real scalar jlo)
{
	real scalar  jm, jhi, inc
	real scalar  ascnd
	real scalar  n

	n = length(xx)
	ascnd=(xx[n] >= xx[1])
	if (jlo <= 0 | jlo > n) {
		jlo = 0
		jhi = n + 1
	}
	else {
		inc = 1
		if (x >= xx[jlo] == ascnd) {
			if (jlo == n) return
			jhi = jlo + 1
			while (x >= xx[jhi] == ascnd) {
				jlo = jhi
				inc = inc + inc
				jhi = jlo + inc
				if (jhi > n) {
					jhi = n + 1
					break
				}
			}
		}
		else {
			if (jlo == 1) {
				jlo = 0
				return
			}
			jhi = jlo--
			while (x < xx[jlo] == ascnd) {
				jhi = jlo
				inc = inc + inc
				if (inc >= jhi) {
					jlo = 0
					break
				}
				else  jlo = jhi - inc
			}
		}
	}
	while (jhi-jlo != 1) {
		jm = trunc((jhi+jlo)/2)
		if (x >= xx[jm] == ascnd)
		 jlo = jm
		else
		 jhi = jm
	}
	if (x == xx[n]) jlo = n - 1
	if (x == xx[1]) jlo = 1
}

end


*! {smcl}
*! {marker mm_clip}{bf:mm_clip.mata}{asis}
*! version 1.0.0  11jul2020  Ben Jann
version 9.2
mata:

real matrix mm_clip(real matrix X, real scalar min, real scalar max, | real scalar miss)
{
    real matrix R
    
    if (args()<4) miss = 0
    if (isfleeting(X)) {
        _mm_clip(X, min, max, miss)
        return(X)
    }
    R = X
    _mm_clip(R, min, max, miss)
    return(R)
}

void _mm_clip(real matrix X, real scalar min, real scalar max, | real scalar miss)
{
    real scalar i, j, r, x
    
    if (args()<4) miss = 0
    if (min>max) _error(3300)
    r = rows(X)
    if (miss) {
        for (j=cols(X); j; j--) {
            for (i=r; i; i--) {
                x = X[i,j]
                if (x>=.) continue
                if (x<min)      X[i,j] = min
                else if (x>max) X[i,j] = max
            }
        }
        return
    }
    for (j=cols(X); j; j--) {
        for (i=r; i; i--) {
            x = X[i,j]
            if (x<min)      X[i,j] = min
            else if (x>max) X[i,j] = max
        }
    }
}

real matrix mm_clipmin(real matrix X, real scalar min)
{
    real matrix R
    
    if (isfleeting(X)) {
        _mm_clipmin(X, min)
        return(X)
    }
    R = X
    _mm_clipmin(R, min)
    return(R)
}

void _mm_clipmin(real matrix X, real scalar min)
{
    real scalar i, j, r
    
    r = rows(X)
    for (j=cols(X); j; j--) {
        for (i=r; i; i--) {
            if (X[i,j]<min) X[i,j] = min
        }
    }
}

real matrix mm_clipmax(real matrix X, real scalar max, | real scalar miss)
{
    real matrix R
    
    if (args()<3) miss = 0
    if (isfleeting(X)) {
        _mm_clipmax(X, max, miss)
        return(X)
    }
    R = X
    _mm_clipmax(R, max, miss)
    return(R)
}

void _mm_clipmax(real matrix X, real scalar max, | real scalar miss)
{
    real scalar i, j, r, x
    
    if (args()<3) miss = 0
    r = rows(X)
    if (miss) {
        for (j=cols(X); j; j--) {
            for (i=r; i; i--) {
                x = X[i,j]
                if (x>=.) continue
                if (x>max) X[i,j] = max
            }
        }
        return
    }
    for (j=cols(X); j; j--) {
        for (i=r; i; i--) {
            if (X[i,j]>max) X[i,j] = max
        }
    }
}

end


*! {smcl}
*! {marker mm_cond}{bf:mm_cond.mata}{asis}
*! version 1.0.0  29feb2008  Ben Jann
version 9.2
mata:

transmorphic matrix mm_cond(real matrix x, transmorphic matrix y, transmorphic matrix z)
{
        transmorphic matrix res
        real scalar r, R, c, C
        transmorphic scalar rx, cx, ry, cy, rz, cz

        if (eltype(y) != eltype(z)) _error(3250)
        R = max((rows(x),rows(y),rows(z)))
        C = max((cols(x),cols(y),cols(z)))
        rx = (rows(x)==1 ? &1 : (rows(x)<R ? _error(3200) : &r))
        cx = (cols(x)==1 ? &1 : (cols(x)<C ? _error(3200) : &c))
        ry = (rows(y)==1 ? &1 : (rows(y)<R ? _error(3200) : &r))
        cy = (cols(y)==1 ? &1 : (cols(y)<C ? _error(3200) : &c))
        rz = (rows(z)==1 ? &1 : (rows(z)<R ? _error(3200) : &r))
        cz = (cols(z)==1 ? &1 : (cols(z)<C ? _error(3200) : &c))
        res = J(R,C, missingof(y))
        for (r=1;r<=R;r++) {
                for (c=1;c<=C;c++) {
                        res[r,c] = (x[*rx,*cx] ? y[*ry,*cy] : z[*rz,*cz])
                }
        }
        return(res)
}

end


*! {smcl}
*! {marker mm_expand}{bf:mm_expand.mata}{asis}
*! version 1.0.0, Ben Jann, 07jun2006
version 9.2
mata:

transmorphic matrix mm_expand(transmorphic matrix X, real vector wr,
 | real vector wc, real scalar sort)
{
	transmorphic matrix Xnew

	if (args()<3) wc = 1
	if (args()<4) sort = 0
	if (isfleeting(X)) {
		_mm_expand(X, wr, wc, sort)
		return(X)
	}
	else {
		_mm_expand(Xnew=X, wr, wc, sort)
		return(Xnew)
	}
}

void _mm_expand(transmorphic matrix X, real vector wr,
 real vector wc, real scalar sort)
{
	real scalar i, b, e, n, add
	real colvector f
	pointer scalar fi

// expand rows
	if (length(wr)==1 & sort==0) _mm_repeat(X, wr, 1)
	else {
		f = trunc(wr:-1) :* (wr:>0)
		_editmissing(f,0)
		n = rows(X)
		if (length(f)==1) {
			fi = &1
			add = n*f
		}
		else {
			if (n!=length(f)) _error(3200)
			fi = &i
			add = sum(f)
		}
		if (add>0) {
			X = X \ J(add,cols(X),missingof(X))
			if (cols(X)>0) {
				if (sort) f = f :+ 1
				b = rows(X)+1
				for (i=n; i>=1; i--) {
					if (f[*fi]<1) continue
					e = b - 1
					b = b - f[*fi]
					X[|b,1 \ e,.|] = X[J(f[*fi],1,i), .]
				}
			}
		}
	}
// expand columns
	if (length(wc)==1 & sort==0) {
		_mm_repeat(X, 1, wc)
		return
	}
	f = trunc(wc:-1) :* (wc:>0)
	_editmissing(f,0)
	n = cols(X)
	if (length(f)==1) {
		fi = &1
		add = n*f
	}
	else {
		if (n!=length(f)) _error(3200)
		fi = &i
		add = sum(f)
	}
	if (add>0) {
		X = X , J(rows(X),add,missingof(X))
		if (rows(X)>0) {
			if (sort) f = f :+ 1
			b = cols(X)+1
			for (i=n; i>=1; i--) {
				if (f[*fi]<1) continue
				e = b - 1
				b = b - f[*fi]
				X[|1,b \ .,e|] = X[., J(1,f[*fi],i)]
			}
		}
	}
}

transmorphic matrix mm_repeat(transmorphic matrix X, real scalar wr,
 | real scalar wc)
{
	transmorphic matrix Xnew

	if (args()<3) wc = 1
	if (isfleeting(X)) {
		_mm_repeat(X, wr, wc)
		return(X)
	}
	else {
		_mm_repeat(Xnew=X, wr, wc)
		return(Xnew)
	}
}

void _mm_repeat(transmorphic matrix X, real scalar wr,
 real scalar wc)
{
	real scalar i, rr, rc, r, c

	r = rows(X)
	c = cols(X)
	rr = max((trunc(wr-1),0))
	if (rr>0) {
		X = X \ J(rr*r, c, missingof(X))
		if (r>0 & c>0)
		 for (i=1; i<=rr; i++) X[|i*r+1,1 \ i*r+r,.|] = X[|1,1 \ r,.|]
	}
	rc = max((trunc(wc-1),0))
	if (rc>0) {
		X = X , J(rows(X), rc*c, missingof(X))
		if (r>0 & c>0)
		 for (i=1; i<=rc; i++) X[|1,i*c+1 \ .,i*c+c|] = X[|1,1 \ .,c|]
	}
}

end


*! {smcl}
*! {marker mm_sort}{bf:mm_sort.mata}{asis}
*! version 1.0.0  09jul2020  Ben Jann
version 9.2
mata:

transmorphic matrix mm_sort(transmorphic matrix X,
    | real rowvector idx, real scalar stable)
{
    if (args()<2) return(X[mm_order(X), .])
    if (args()<3) return(X[mm_order(X, idx), .])
    return(X[mm_order(X, idx, stable), .])
}

real colvector mm_order(transmorphic matrix X,
    | real rowvector idx, real scalar stable)
{
    real scalar r, c, l
    
    if (args()<3) stable = 0
    // redirect to order()
    c = cols(X)
    if (!stable) {
        if (args()>=2 & idx!=.) return(order(X, idx))
        return(order(X, (c>0 ? 1..c : J(1,0,.))))
    }
    // check idx
    l = (idx==. ? 0 : length(idx))
    if (l>c) _error(3200)
    // one row or less: nothing to do
    r = rows(X)
    if (r<=1) return(J(r,1,1))
    // zero columns: return index
    if (c==0) return(1::r)
    // X is real: append index for stable order
    if (isreal(X)) {
        if (l==0) return(order((X,(1::r)), 1..c+1))
        if (anyof(idx,0)) _error(3300)
        return(order((X[,abs(idx)],(1::r)), ((1..l):*sign(idx), l+1)))
    }
    // X is not real; need to create a group id
    if (l==0) return(order((mm_group(X), (1::r)), (1,2)))
    return(order((mm_group(X, idx), (1::r)), (1,2)))
}

end



*! {smcl}
*! {marker mm_unorder2}{bf:mm_unorder2.mata}{asis}
*! version 1.0.0, Ben Jann, 29mar2006
version 9.0
mata:

real colvector mm_unorder2(real scalar n) return(order(uniform(n,2),(1,2)))

end


*! {smcl}
*! {marker mm_jumble2}{bf:mm_jumble2.mata}{asis}
*! version 1.0.0, Ben Jann, 29mar2006
version 9.0
mata:

transmorphic matrix mm_jumble2(transmorphic matrix x)
{
	return(x[mm_unorder2(rows(x)), .])
}

end


*! {smcl}
*! {marker mm__jumble2}{bf:mm__jumble2.mata}{asis}
*! version 1.0.0, Ben Jann, 29mar2006
version 9.0
mata:

function mm__jumble2(transmorphic matrix x)
{
	_collate(x,mm_unorder2(rows(x)))
}

end


*! {smcl}
*! {marker mm_pieces}{bf:mm_pieces.mata}{asis}
*! version 1.0.2, Ben Jann, 01jun2015
version 9.2
mata:

string rowvector mm_pieces(
    string scalar s,
    real scalar w,
    | real scalar nobreak,
      real scalar ascii)
{
    real scalar         i, j, k, n, l, b, nobr
    string scalar c
    string rowvector    res

    nobr = ( args()>2 ? nobreak : 0 )
    if (stataversion()>=1400 & !(args()>3 ? ascii : 0)) {
        return(_mm_pieces_goto14(s, w, nobr))
    }
    l = strlen(s)
    if (l<2 | w>=l) return(strtrim(s))
    res = J(1, _mm_npieces(s, w, nobr, (args()>3 ? ascii : 0)), "")
    j = k = n = 0
    b = 1
    for (i=1; i<=l; i++) {
        c = substr(s, i, 1)
        if (j<1) { // skip to first nonblank character
            if (c==" ") {
                b++
                continue
            }
        }
        j++
        if (i==l) res[++n] = strtrim(substr(s, b, .))
        else {
            if (c==" ") k = i
            if (j>=w) {
                if (k<1) {
                    if (nobr) continue
                    k = i
                }
                else {
                    if (substr(s, i+1, 1)==" ") k = i
                }
                res[++n] = strtrim(substr(s, b, k-b+1))
                j = i - k
                b = k + 1
                k = 0
            }
        }
    }
    return(res)
}

real matrix mm_npieces(
    string matrix S,
    real scalar w,
    | real scalar nobreak,
      real scalar ascii)
{
    real scalar  i, j
    real matrix  res

    res = J(rows(S),cols(S),1)
    if (w>=.) return(res)
    for (i=1; i<=rows(S); i++) {
        for (j=1; j<=cols(S); j++) {
            res[i,j] = _mm_npieces(S[i,j], w, (args()>2 ? nobreak : 0), 
                                              (args()>3 ? ascii : 0))
        }
    }
    return(res)
}

real scalar _mm_npieces(
    string scalar s,
    real scalar w,
    | real scalar nobreak,
      real scalar ascii)
{
    real scalar   i, j, k, n, l, nobr
    string scalar c

    nobr = ( args()>2 ? nobreak : 0 )
    if (stataversion()>=1400 & !(args()>3 ? ascii : 0)) {
        return(_mm_npieces_goto14(s, w, nobr))
    }
    l = strlen(s)
    if (l<2 | w>=l) return(1)
    j = k = n = 0
    for (i=1; i<=l; i++) {
        c = substr(s, i, 1)
        if (j<1) { // skip to first nonblank character
            if (c==" ") continue
        }
        j++
        if (i==l) n++
        else {
            if (c==" ") k = i
            if (j>=w) {
                if (k<1) {
                    if (nobr) continue
                    k = i
                }
                else {
                    if (substr(s, i+1, 1)==" ") k = i
                }
                n++
                j = i - k
                k = 0
            }
        }
    }
    if (n==0) n++
    return(n)
}

string rowvector _mm_pieces_goto14(string scalar s,
    real scalar w,
    real scalar nobr)
{
    return(_mm_pieces14(s, w, nobr))
}

real scalar _mm_npieces_goto14(string scalar s,
    real scalar w,
    real scalar nobr)
{
    return(_mm_npieces14(s, w, nobr))
}

end


*! {smcl}
*! {marker mm_regexr}{bf:mm_regexr.mata}{asis}
*! version 1.0.0  31jan2017  Ben Jann
version 9.2
mata:

string matrix mm_regexr(string matrix x, string matrix y, string matrix z)
{
    string matrix       res
    real scalar         r, R, c, C
    transmorphic scalar rx, cx, ry, cy, rz, cz
    pragma unset        r
    pragma unset        c

    R = max((rows(x),rows(y),rows(z)))
    C = max((cols(x),cols(y),cols(z)))
    rx = (rows(x)==1 ? &1 : (rows(x)<R ? _error(3200) : &r))
    cx = (cols(x)==1 ? &1 : (cols(x)<C ? _error(3200) : &c))
    ry = (rows(y)==1 ? &1 : (rows(y)<R ? _error(3200) : &r))
    cy = (cols(y)==1 ? &1 : (cols(y)<C ? _error(3200) : &c))
    rz = (rows(z)==1 ? &1 : (rows(z)<R ? _error(3200) : &r))
    cz = (cols(z)==1 ? &1 : (cols(z)<C ? _error(3200) : &c))
    res = J(R,C, "")
    for (r=1;r<=R;r++) {
        for (c=1;c<=C;c++) {
            res[r,c] = _mm_regexr(x[*rx,*cx], y[*ry,*cy], z[*rz,*cz])
        }
    }
    return(res)
}

string scalar _mm_regexr(string scalar s0, string scalar from, string scalar to)
{
    real scalar         i, j
    string scalar       s, t, BSLASH
    string rowvector    sub
    pragma unset        s
    
    if (regexm(s0, from)) {
        sub = regexs()
        sub = (sub, J(1, 10-cols(sub), ""))
        BSLASH = "\"
        for (i=1; i<=strlen(to); i++) {
            if ((t=substr(to,i,1))==BSLASH) {
                i++
                if ((t=substr(to,i,1))==BSLASH) s = s + BSLASH      // "\\"
                else if ((j=strtoreal(t))<.)    s = s + sub[j+1]    // "\#"
                else                            s = s + BSLASH + t
            }
            else s = s + t
        }
        return(subinstr(s0, sub[1], s, 1))
    }
    else return(s0)
}

end


*! {smcl}
*! {marker mm_invtokens}{bf:mm_invtokens.mata}{asis}
*! version 1.0.2  Ben Jann  05jan2008
version 9.0
mata:

string scalar mm_invtokens(string vector In, | real scalar noclean)
{
    string scalar Out
    string scalar tok
    real scalar i

    if (args()<2) noclean = 0
    Out = ""
    for (i=1; i<=length(In); i++) {
        tok = In[i]
        if (noclean)                 tok = "`" + `"""' + tok + `"""' + "'"
        else if (strpos(tok, `"""')) tok = "`" + `"""' + tok + `"""' + "'"
        else if (strpos(tok, " "))   tok = `"""' + tok + `"""'
        else if (tok=="")            tok = `"""' + `"""'
        if (i>1) tok = " " + tok
        Out = Out + tok
    }
    return(Out)
}
end


*! {smcl}
*! {marker mm_realofstr}{bf:mm_realofstr.mata}{asis}
*! version 1.0.2, Ben Jann, 24mar2006
version 9.0
mata:

real matrix mm_realofstr(string matrix S)
{
	real matrix R
	string scalar tmp
	real scalar i, j

	R = J(rows(S), cols(S), .)
	tmp = st_tempname()
	for (i=1; i<=rows(S); i++) {
		for (j=1; j<=cols(S); j++) {
			stata("scalar " + tmp + "=real(`" + `"""' + S[i,j] + `"""' + "')")
			R[i,j] = st_numscalar(tmp)
		}
	}
	stata("capture scalar drop " + tmp)
	return(R)
}
end


*! {smcl}
*! {marker mm_strexpand}{bf:mm_strexpand.mata}{asis}
*! version 1.0.4, Ben Jann, 30apr2007
version 9.0
mata:

string scalar mm_strexpand(string scalar s, string vector slist,
 | string scalar def, real scalar unique, string scalar errtxt)
{
	real scalar   err
	string scalar res

	if (args()<5) errtxt = `"""' + s + `"" invalid"'
	if (args()<4) unique = 0
	err = _mm_strexpand(res, s, slist, def, unique)
	if (err) _error(err, errtxt)
	return(res)
}

real scalar _mm_strexpand(res, string scalar s,
 string vector slist, | string scalar def, real scalar unique)
{
	real scalar i, l, match

	if (s=="") {
		res = def
		return(0)
	}
	if (args()<5) unique = 0
	l = strlen(s)
	if (unique) {
		match = 0
		for (i=1; i<=length(slist); i++) {
			if (s==substr(slist[i], 1, l)) {
				if (match) return(3498)
				match = i
			}
		}
		if (match) {
			res = slist[match]
			return(0)
		}
	}
	else {
		for (i=1; i<=length(slist); i++) {
			if (s==substr(slist[i], 1, l)) {
				res = slist[i]
				return(0)
			}
		}
	}
	return(3499)
}
end


*! {smcl}
*! {marker mm_matlist}{bf:mm_matlist.mata}{asis}
*! version 1.0.4  31aug2007  Ben Jann
version 9.2
mata:

void mm_matlist(
    real matrix X,
    | string matrix fmt,
      real scalar b,
      string vector rstripe,
      string vector cstripe,
      string scalar rtitle,
      transmorphic matrix res
    )
{
    real scalar             i, j, rw, lw, j0, j1, l, r, save
    real rowvector          wd
    string rowvector        sfmt
    pointer(real scalar) scalar     fi, fj
    pointer(string vector) scalar   rstr, cstr

    if (args()<2) fmt = "%9.0g"
    if (args()<3) b = 1
    save = (args()==7)
    //wd = strlen(sprintf(fmt,-1/3))

    if ((rows(fmt)!=1 & rows(fmt)!=rows(X)) |
        (cols(fmt)!=1 & cols(fmt)!=cols(X))) _error(3200)
    if (!anyof((0,1,2,3),b)) _error(3300)

    fi = (rows(fmt)!=1 ? &i : &1)
    fj = (cols(fmt)!=1 ? &j : &1)

    wd = J(1,cols(X),0)
    for (i=1;i<=rows(X);i++) {
        for (j=1;j<=cols(X);j++) {
            wd[j] = max((wd[j],strlen(sprintf(fmt[*fi,*fj],X[i,j]))))
        }
    }

    if (length(X)==0) return

    if (length(rstripe)<1)  rstr = &strofreal(1::rows(X))
    else rstr = &rstripe
    if (length(cstripe)<1)  cstr = &strofreal(1..cols(X))
    else if (rows(cstripe)!=1)  cstr = &(cstripe')
    else                        cstr = &cstripe

    if ((*rstr!="" & length(*rstr)!=rows(X)) |
        (*cstr!="" & length(*cstr)!=cols(X))) _error(3200)

    if (*rstr!="")  rw = max((max(strlen(*rstr)),strlen(rtitle)))
    else            rw = -1
    if (*cstr!="")  wd = colmax((strlen(*cstr)\wd)) :+ 2
    else            wd = wd :+ 2

    sfmt = "%" :+ strofreal(wd) :+ "s"

    if (save) {
        _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle,
            wd, rw, sfmt, 1, cols(X), 1, 1, res)
        return
    }

    lw = st_numscalar("c(linesize)") - ((2+rw+2*(b>0))+2*(b==1))
    if ((sum(wd)+cols(X))<=lw) {
        _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle,
            wd, rw, sfmt, 1, cols(X), 1, 1)
        return
    }

    j0 = 1
    l = 1
    r = 0
    while (j0<=cols(X)) {
        j1 = j0
        while (sum(wd[|j0\j1|]:+1)<=lw) {
            if (++j1>cols(X)) {
                r = 1
                break
            }
        }
        j1 = max((j1-1,j0))
        _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle,
            wd, rw, sfmt, j0, j1, l, r)
        l = 0
        j0 = j1+1
    }
}

string colvector mm_smatlist(
    real matrix X,
    | string matrix fmt,
      real scalar b,
      string vector rstripe,
      string vector cstripe,
      string scalar rtitle
    )
{
    string colvector res

    if (args()<2) fmt = "%9.0g"
    if (args()<3) b = 1
    if (args()<4) rstripe = J(0,1,"")
    if (args()<5) cstripe = J(1,0,"")
    if (args()<5) rtitle = ""
    mm_matlist(X, fmt, b, rstripe, cstripe, rtitle, res)
    return(res)
}

void _mm_matlist(
    real matrix X,
    string matrix fmt,
    real scalar b,
    string vector rstripe,
    string vector cstripe,
    string scalar rtitle,
    real rowvector wd,
    real scalar rw,
    string rowvector sfmt,
    real scalar j0,
    real scalar j1,
    real scalar l,
    real scalar r,
    | transmorphic matrix res
    )
{
    real scalar          i, j, di, ri
    string scalar        rfmt, line
    pointer(real scalar) scalar fi, fj

    di = args()<14
    fi = (rows(fmt)!=1 ? &i : &1)
    fj = (cols(fmt)!=1 ? &j : &1)
    rfmt = "%"+strofreal(rw)+"s"

    if (di==0) {
        ri = 0
        res = J(rows(X)+(1+(b==3))*(cstripe!="")+(b>0)+(b==1)
            +(b==3)+(l==0&(b==0|b==2)),1,"")
    }

    if (l==0 & (b==0 | b==2)) {
        if (di) printf("\n")
        else res[++ri] = ""
    }
    if (cstripe!="") {
        if (b==3) {
            line = sprintf("{txt}" + "{hline " + strofreal(2+rw+1) +  "}{c TT}"
                + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)-1) + "}")
            if (di) display(line)
            else res[++ri] = line
        }
        line = sprintf("{txt}  ")
        if (rstripe!="") line = line + sprintf(rfmt+" ", rtitle)
        if (b>0) line = line + sprintf((b==1 ? " " : "{c |}"))
        for (j=j0;j<=j1;j++) {
            line = line + sprintf(sfmt[j], cstripe[j])
            if (j<j1) line = line + sprintf(" ")
        }
        if (di) display(line)
        else res[++ri] = line
    }
    if (b) {
        line = sprintf("{txt}" + (b==1 ? (2+rw+1)*" " +
            (l ? "{c TLC}" : " ") : "{hline " + strofreal(2+rw+1) + "}{c +}")
            + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)+1-2*(b!=1)) + "}" +
            (b==1 ? (r ? "{c TRC}" : "") : ""))
        if (di) display(line)
        else res[++ri] = line
    }
    for (i=1;i<=rows(X);i++) {
        line = sprintf("{txt}  ")
        if (rstripe!="") line = line + sprintf(rfmt+" ", rstripe[i])
        if (b>0) line = line + sprintf((l | b!=1 ? "{c |}" : " "))
        line = line + sprintf("{res}")
        for (j=j0;j<=j1;j++) {
            line = line + sprintf(sfmt[j],sprintf(fmt[*fi,*fj], X[i,j]))
            if (j<j1) line = line + sprintf(" ")
        }
        if (b==1 & r) line = line + sprintf("  {txt}{c |}")
        if (di) display(line)
        else res[++ri] = line
    }
    if (b==1|b==3) {
        line = sprintf("{txt}" + (b==1 ? (2+rw+1)*" " +
            (l ? "{c BLC}" : " ") : "{hline " + strofreal(2+rw+1) + "}{c BT}")
            + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)+1-2*(b!=1)) + "}" +
            (b==1 ? (r ? "{c BRC}" : "") : ""))
        if (di) display(line)
        else res[++ri] = line
    }
}

end


*! {smcl}
*! {marker mm_insheet}{bf:mm_insheet.mata}{asis}
*! version 2.0.1, Ben Jann, 18mar2006
*! based on cat(), version 2.0.0  23jan2006
version 9.0
mata:

string matrix mm_insheet(string scalar filename,
 | string scalar del, real scalar uline1, real scalar line2)
{
	real scalar    fh, fpos, line1, i, r, c, j, delpos, ji
	string matrix  EOF
	string matrix  res
	string scalar  line

	/* ------------------------------------------------------------ */
	/* setup */
	if (del=="") del = char(9)
	fh  = fopen(filename, "r")
	EOF = J(0, 0, "")
	line1 = floor(uline1)
	if (line1<1 | line1>=.) line1 = 1

	/* ------------------------------------------------------------ */
	/* nothing to read case  */
	if (line1>line2) {
		fclose(fh)
		return(J(0, 0, ""))
	}

	/* ------------------------------------------------------------ */
	 /* skip opening lines  */
	for (i=1; i<line1; i++) {
		if (fget(fh)==EOF) {
			fclose(fh)
			return(J(0, 0, ""))
		}
	}
	fpos = ftell(fh)

	/* ------------------------------------------------------------ */
	/* count lines and columns  */
	j = 0
	for (i=line1; i<=line2; i++) {
		if ((line=fget(fh))==EOF) break
		ji = 1
		while (delpos = strpos(line, del)) {
			ji++
			line = substr(line, delpos+1, .)
		}
		if (ji>j) j = ji
	}
	res = J(r = i-line1, c = j, "")

	/* ------------------------------------------------------------ */
	/* read lines */
	fseek(fh, fpos, -1)
	for (i=1; i<=r; i++) {
		if ((line=fget(fh))==EOF) {
			/* unexpected EOF -- file must have changed */
			fclose(fh)
			if (--i) return(res[|1,1 \ i,c|])
			return(J(0,c,""))
		}
		for (j=1; j<=c; j++) {
			delpos = strpos(line, del)
			if (delpos==0) {
				res[i,j] = line
				break
			}
			res[i,j] = substr(line, 1, delpos-1)
			line = substr(line, delpos+1, .)
		}
	}
	fclose(fh)
	return(res)
}

end


*! {smcl}
*! {marker mm_infile}{bf:mm_infile.mata}{asis}
*! version 2.0.2, Ben Jann, 18mar2006
*! based on cat(), version 2.0.0  23jan2006
version 9.0
mata:

string matrix mm_infile(string scalar filename,
 | real scalar uline1, real scalar line2)
{
	real scalar    fh, fpos, line1, i, r, c, j
	string matrix  EOF
	string matrix  res
	string scalar  line

	/* ------------------------------------------------------------ */
	/* setup */
	fh  = fopen(filename, "r")
	EOF = J(0, 0, "")
	line1 = floor(uline1)
	if (line1<1 | line1>=.) line1 = 1

	/* ------------------------------------------------------------ */
	/* nothing to read case  */
	if (line1>line2) {
		fclose(fh)
		return(J(0, 0, ""))
	}

	/* ------------------------------------------------------------ */
	 /* skip opening lines  */
	for (i=1; i<line1; i++) {
		if (fget(fh)==EOF) {
			fclose(fh)
			return(J(0, 0, ""))
		}
	}
	fpos = ftell(fh)

	/* ------------------------------------------------------------ */
	/* count lines and columns  */
	j = 0
	for (i=line1; i<=line2; i++) {
		if ((line=fget(fh))==EOF) break
		line = tokens(line)
		if (cols(line)>j) j = cols(line)
	}
	res = J(r = i-line1, c = j, "")

	/* ------------------------------------------------------------ */
	/* read lines */
	fseek(fh, fpos, -1)
	for (i=1; i<=r; i++) {
		if ((line=fget(fh))==EOF) {
			/* unexpected EOF -- file must have changed */
			fclose(fh)
			if (--i) return(res[|1,1 \ i,c|])
			return(J(0,c,""))
		}
		line = tokens(line)
		res[|i,1 \ i,cols(line)|] = line
	}
	fclose(fh)
	return(res)
}

end


*! {smcl}
*! {marker mm_outsheet}{bf:mm_outsheet.mata}{asis}
*! version 1.0.4, Ben Jann, 14apr2006
version 9.0
mata:

void mm_outsheet(string scalar fn, string matrix s,
 | string scalar mode0, string scalar del)
{
	string scalar line, mode, m
	real scalar i, j, fh

	if (args()<4) del = char(9)
	mode = mm_strexpand(mode0, ("append", "replace"))
	m = "w"
	if (mode=="replace") unlink(fn)
	else if (mode=="append") m = "a"
	fh = fopen(fn, m)
	for (i=1; i<=rows(s); i++) {
		line = J(1,1,"")
		for (j=1; j<=cols(s); j++) {
			line = line + s[i,j]
			if (j<cols(s)) line = line + del
		}
		fput(fh, line)
	}
	fclose(fh)
}
end


*! {smcl}
*! {marker mm_callf}{bf:mm_callf.mata}{asis}
*! version 1.0.1, Ben Jann, 16jun2006
version 9.2

local o1 "o1"
local po1 "*p.o1"
forv i=2/10 {
	local o`i' "`o`=`i'-1'', o`i'"
	local po`i' "`po`=`i'-1'', *p.o`i'"
}

mata:

struct mm_callf_o10 {
	pointer(function) scalar f
	real scalar              n
	pointer scalar           `o10'
}

struct mm_callf_o10 scalar mm_callf_setup(
 pointer(function) scalar f, real scalar n, | `o10')
{
	struct mm_callf_o10 scalar p

	p.f  = f
	p.n  = n
	p.o1 = &o1 ; p.o2 = &o2 ; p.o3 = &o3 ; p.o4 = &o4 ; p.o5 = &o5
	p.o6 = &o6 ; p.o7 = &o7 ; p.o8 = &o8 ; p.o9 = &o9 ; p.o10 = &o10
	return(p)
}

transmorphic mm_callf(struct mm_callf_o10 p, | `o5')
{
	if (args()==1) return(_mm_callf0(p))
	if (args()==2) return(_mm_callf1(p, `o1'))
	if (args()==3) return(_mm_callf2(p, `o2'))
	if (args()==4) return(_mm_callf3(p, `o3'))
	if (args()==5) return(_mm_callf4(p, `o4'))
	               return(_mm_callf5(p, `o5'))
}

transmorphic _mm_callf0(struct mm_callf_o10 p)
{
	if (p.n<=0) return( (*p.f)() )
	if (p.n==1) return( (*p.f)(`po1') )
	if (p.n==2) return( (*p.f)(`po2') )
	if (p.n==3) return( (*p.f)(`po3') )
	if (p.n==4) return( (*p.f)(`po4') )
	if (p.n==5) return( (*p.f)(`po5') )
	if (p.n==6) return( (*p.f)(`po6') )
	if (p.n==7) return( (*p.f)(`po7') )
	if (p.n==8) return( (*p.f)(`po8') )
	if (p.n==9) return( (*p.f)(`po9') )
	            return( (*p.f)(`po10') )
}

transmorphic _mm_callf1(struct mm_callf_o10 p, `o1')
{
	if (p.n<=0) return( (*p.f)(`o1') )
	if (p.n==1) return( (*p.f)(`o1', `po1') )
	if (p.n==2) return( (*p.f)(`o1', `po2') )
	if (p.n==3) return( (*p.f)(`o1', `po3') )
	if (p.n==4) return( (*p.f)(`o1', `po4') )
	if (p.n==5) return( (*p.f)(`o1', `po5') )
	if (p.n==6) return( (*p.f)(`o1', `po6') )
	if (p.n==7) return( (*p.f)(`o1', `po7') )
	if (p.n==8) return( (*p.f)(`o1', `po8') )
	if (p.n==9) return( (*p.f)(`o1', `po9') )
	            return( (*p.f)(`o1', `po10') )
}

transmorphic _mm_callf2(struct mm_callf_o10 p, `o2')
{
	if (p.n<=0) return( (*p.f)(`o2') )
	if (p.n==1) return( (*p.f)(`o2', `po1') )
	if (p.n==2) return( (*p.f)(`o2', `po2') )
	if (p.n==3) return( (*p.f)(`o2', `po3') )
	if (p.n==4) return( (*p.f)(`o2', `po4') )
	if (p.n==5) return( (*p.f)(`o2', `po5') )
	if (p.n==6) return( (*p.f)(`o2', `po6') )
	if (p.n==7) return( (*p.f)(`o2', `po7') )
	if (p.n==8) return( (*p.f)(`o2', `po8') )
	if (p.n==9) return( (*p.f)(`o2', `po9') )
	            return( (*p.f)(`o2', `po10') )
}

transmorphic _mm_callf3(struct mm_callf_o10 p, `o3')
{
	if (p.n<=0) return( (*p.f)(`o3') )
	if (p.n==1) return( (*p.f)(`o3', `po1') )
	if (p.n==2) return( (*p.f)(`o3', `po2') )
	if (p.n==3) return( (*p.f)(`o3', `po3') )
	if (p.n==4) return( (*p.f)(`o3', `po4') )
	if (p.n==5) return( (*p.f)(`o3', `po5') )
	if (p.n==6) return( (*p.f)(`o3', `po6') )
	if (p.n==7) return( (*p.f)(`o3', `po7') )
	if (p.n==8) return( (*p.f)(`o3', `po8') )
	if (p.n==9) return( (*p.f)(`o3', `po9') )
	            return( (*p.f)(`o3', `po10') )
}

transmorphic _mm_callf4(struct mm_callf_o10 p, `o4')
{
	if (p.n<=0) return( (*p.f)(`o4') )
	if (p.n==1) return( (*p.f)(`o4', `po1') )
	if (p.n==2) return( (*p.f)(`o4', `po2') )
	if (p.n==3) return( (*p.f)(`o4', `po3') )
	if (p.n==4) return( (*p.f)(`o4', `po4') )
	if (p.n==5) return( (*p.f)(`o4', `po5') )
	if (p.n==6) return( (*p.f)(`o4', `po6') )
	if (p.n==7) return( (*p.f)(`o4', `po7') )
	if (p.n==8) return( (*p.f)(`o4', `po8') )
	if (p.n==9) return( (*p.f)(`o4', `po9') )
	            return( (*p.f)(`o4', `po10') )
}

transmorphic _mm_callf5(struct mm_callf_o10 p, `o5')
{
	if (p.n<=0) return( (*p.f)(`o5') )
	if (p.n==1) return( (*p.f)(`o5', `po1') )
	if (p.n==2) return( (*p.f)(`o5', `po2') )
	if (p.n==3) return( (*p.f)(`o5', `po3') )
	if (p.n==4) return( (*p.f)(`o5', `po4') )
	if (p.n==5) return( (*p.f)(`o5', `po5') )
	if (p.n==6) return( (*p.f)(`o5', `po6') )
	if (p.n==7) return( (*p.f)(`o5', `po7') )
	if (p.n==8) return( (*p.f)(`o5', `po8') )
	if (p.n==9) return( (*p.f)(`o5', `po9') )
	            return( (*p.f)(`o5', `po10') )
}

end
